# Ion Clusters and Networks in Water-in-Salt Electrolytes

Michael McEldrew  Department of Chemical Engineering, Massachusetts Institute
of Technology, Cambridge, MA, USA  Zachary A. H. Goodwin  Department of
Physics, CDT Theory and Simulation of Materials, Imperial College of London,
South Kensington Campus, London SW7 2AZ, UK  Department of Chemistry, Imperial
College of London, South Kensington Campus, London SW7 2AZ, UK  Sheng Bi
State Key Laboratory of Coal Combustion, School of Energy and Power
Engineering, Huazhong University of Science and Technology (HUST), Wuhan,
Hubei 430074, China  Alexei A. Kornyshev  Department of Chemistry, Imperial
College of London, Molecular Sciences Research Hub, White City Campus, Wood
Lane, London W12 0BZ, UK  Thomas Young Centre for Theory and Simulation of
Materials, Imperial College of London, South Kensington Campus, London SW7
2AZ, UK  Institute of Molecular Science and Engineering, Imperial College of
London, South Kensington Campus, London SW7 2AZ, UK  Martin Z. Bazant
<EMAIL_ADDRESS>Department of Chemical Engineering, Massachusetts Institute of
Technology, Cambridge, MA, USA  Department of Mathematics, Massachusetts
Institute of Technology, Cambridge, MA, USA

###### Abstract

Water-in-salt electrolytes (WiSEs) are a class of super-concentrated
electrolytes that have shown much promise in replacing organic electrolytes in
lithium-ion batteries. At the extremely high salt concentrations of WiSEs,
ionic association is more complicated than the simple ion pair description. In
fact, large branched clusters can be present in WiSEs, and past a critical
salt concentration, an infinite percolating ionic network can form
spontaneously. In this work, we simplify our recently developed thermodynamic
model of reversible ionic aggregation and gelation, tailoring it specifically
for WiSEs. Our simplified theory only has a handful of parameters, all of
which can be readily determined from simulations. Our model is able to
quantitatively reproduce the populations of ionic clusters of different sizes
as a function of salt concentration, the critical salt concentration for ionic
gelation, and the fraction of ions incorporated into the ionic gel, as
observed from molecular simulations of three different lithium-based WiSEs.
The extent of ionic association and gelation greatly affects the effective
ionic strength of solution, the coordination environment of active cations
that is known to govern the chemistry of the solid-electrolyte interface, and
the thermodynamic activity of all species in the electrolyte.

##  I Introduction

Recently developed Water-in-Salt Electrolytes (WiSEs) have received a great
deal of attention due to their potential for application in lithium-ion
batteries (LIBs) Suo  _et al._ (2013); Sodeyama  _et al._ (2014); Suo  _et
al._ (2015); Yamada  _et al._ (2016); Wang  _et al._ (2016); Gambou-Bosca and
Bélanger (2016); Sun  _et al._ (2017); Suo  _et al._ (2017); Dong  _et al._
(2017a); Diederichsen  _et al._ (2017); Yang  _et al._ (2017a, b); Wang  _et
al._ (2018); Leonard  _et al._ (2018); Yamada  _et al._ (2019); Yang  _et al._
(2019); Dou  _et al._ (2019); Lewis  _et al._ (2020). WiSEs are super-
concentrated electrolytes often containing a mere 2-3 water molecules per ion
pair. Typically, WiSEs contain small, alkali cations (such as Li+) and bulky,
fluorinated anions [such as TFSI-, Bis(trifluoromethane)sulfonimide]. Owing to
their high concentration, WiSEs have displayed enhanced electrochemical
stability windows (ESW) compared to their dilute counterparts, making them
promising candidates to replace the flammable organic electrolytes that are
widely used in LIBs. Additionally, WiSEs tend to have superior transport
properties compared to organic LIB electrolytes Borodin  _et al._ (2017a). We
direct the reader to refs. 20; 21; 22 for reviews highlighting recent progress
in the literature on WiSEs.

There have been a number of computational studies highlighting the interfacial
Vatamanu and Borodin (2017); McEldrew  _et al._ (2018); Li  _et al._ (2018)
and bulk Borodin  _et al._ (2017a); Choi  _et al._ (2018a); Lim  _et al._
(2018); Yu  _et al._ (2020); Lewis  _et al._ (2020); Andersson  _et al._
(2020) behavior of WiSEs, which have given valuable insight into the
physicochemical and electrochemical behavior of WiSEs. However, WiSEs have not
received enough theoretical attention in the literature. From a theoretical
point of view, WiSE’s embody the middle ground between solvent-less ionic
liquids (ILs) and traditional solvent-dominated electrolytes. Theories of
traditional electrolytes are well-established and have largely been rooted in
the classical theories of Debye and Hückel Debye and Hückel (1923); Hückel
(1925). In the past decade, ILs have received increasing attention from a
theoretical point of view as highly correlated ionic matter as well as variety
of their applications. Fedorov and Kornyshev (2014); Dong  _et al._ (2017b).
However, theories of traditional electrolytes and ILs often require simplified
pictures of ion solvation and ionic aggregation. For example, in traditional
electrolytes the solvent is often modeled as a dielectric continuum that is
dependent on salt concentration and/or electrostatic fields Hückel (1925);
Fawcett and Tikanen (1996); Abbas  _et al._ (2007); Vincze  _et al._ (2010);
Shilov and Lyashchenko (2015), and ion association has been largely limited to
the formation of ion pairs Bjerrum (1926); Ebeling and Grigo (1980); Levin and
Fisher (1996); Marcus and Hefter (2006). Prior to Ref. 42, ionic aggregation
in ILs was treated only through the remaining amount of “free” ions not in
ionic aggregates Chen  _et al._ (2018); Feng  _et al._ (2019); Goodwin and
Kornyshev (2017); Goodwin  _et al._ (2017) or through the formation of ion
pairs Ma  _et al._ (2015); Lee  _et al._ (2014), and in rare cases small
clusters of ion pairs Avni  _et al._ (2020). For WiSE’s, the solvent is
neither abundant nor negligible, but rather it is a discrete yet essential
component of the solution that should be considered on the same footing as
ions in the electrolyte. Similarly, the coordinating environment an ion is
neither a pristine solvation shell of solvent (dilute electrolyte) nor
counter-ion (ionic liquid/crystal), but rather a complex combination of
counter-ion and solvent coordination Pivnic  _et al._ (2019). Moreover, ions
can be clustered in complex, branched aggregates that may even percolate
throughout the mixture Choi  _et al._ (2017); Borodin  _et al._ (2017a); Choi
_et al._ (2018a); Lim  _et al._ (2018); Lewis  _et al._ (2020); McEldrew  _et
al._ (2020b).

In this work, we apply our recently developed theory of ionic aggregation and
gelation of super-concentrate electrolytes McEldrew  _et al._ (2020b) to model
WiSEs. The model accounts for the competition between ion association and ion
solvation, finite molecular volumes, and the free energy to form ionic
associations. We use our model to analyze ion aggregation, solvation, and
gelation in various lithium-based aqueous electrolytes for salt concentrations
ranging from the dilute regime to the water-in-salt regime. Using molecular
dynamics (MD) simulations, we calculate the extent of ionic aggregation and
solvation, distribution of ionic cluster sizes, and even the presence and
extent of percolating ion networks. Our theory only depends on a handful of
physically transparent parameters, such as the volume of ions and the number
of “bonds” that they can form and their strength. The performed MD simulations
were also used to calculate these parameters. Overall, we find excellent
agreement between our theory (with parameters obtained from MD) and the
molecular simulations, which gives one confidence in the predictive power of
the theory.

##  II Theory

In our previous works, we developed a general theory of ion aggregation in
super-concentrated electrolytes McEldrew  _et al._ (2020b), and a particular
application of which was to ILs McEldrew  _et al._ (2020a). The general theory
McEldrew  _et al._ (2020b) accounted for associations between ions, forming
complex and potentially even percolating ion networks, as well as the
hydration of ions by solvent molecules, with bonds between solvent molecules
being neglected.

Here, we simplify that general theory McEldrew  _et al._ (2020b) by making
several assumptions, which will be largely valid for super-concentrated WiSEs:

  1. 1.

Selective bonding of water. We neglect anion solvation, as the anions used in
WiSE’s are bulky and relatively hydrophobic. Since the bonds typically form to
the oxygens of the anions (as shown later), a consistent approximation is to
also neglect bonding between water molecules. The water is only able to bind
to the lithium cations, which means that secondary hydration shells are not
accounted for. While this is certainly a simplified picture of the solution
structure, we believe the Li+-OH2 coordination to be the dominant effect.

  2. 2.

Simplified electrostatics and interactions. We neglect the excess
electrostatic energy of the mixture, such as the Debye-Hückle term of dilute
electrolytes Debye and Hückel (1923); Hückel (1925). This assumption relies on
the electrostatic energy being dominated by the formation of ionic clusters,
which are assumed to behave ideally in the solution (no inter-cluster
enthalpic interactions). This simplification was also employed successfully in
Ref. 42 to model ILs, and we expect it to work reasonably well for super-
concentrated WiSEs. Moreover, in Ref. 52, the effect of Debye-Hückle screening
and Born solvation of free ions was studied and was found to not significantly
affect ion cluster distributions at high concentrations.

  3. 3.

“Sticky” cation approximation. Our last simplification is that the cations
always have a full coordination shell, comprising of water molecules and/or
anions. As we will show, this assumption is especially valid for lithium-based
WiSEs, in which lithium cations interact very strongly with both the water and
the anion.

For brevity, we will simply outline the crux of theory in the main text, and
relegate some of the details to the supporting information (SI).

We consider a polydisperse mixture of $\sum_{lms}N_{lms}$ ionic clusters
containing $l$ cations, $m$ anions, and $s$ water molecules (a rank $lms$
cluster), and an interpenetrating gel network (if present) containing
$N_{+}^{gel}$ cations, $N_{-}^{gel}$ anions, and $N_{0}^{gel}$ water
molecules. As in Refs. 52 & 42, we use a Flory-like lattice fluid free energy,
which has been used extensively to describe associative polymer mixtures
Tanaka (1989, 1990); Tanaka and Stockmayer (1994); Tanaka and Ishida (1995);
Ishida and Tanaka (1997); Tanaka (1998); Tanaka and Ishida (1999); Tanaka
(2002):

$\displaystyle\beta\Delta F$
$\displaystyle=\sum_{l,m,s}\left[N_{lms}\log\left(\phi_{lms}\right)+N_{lms}\Delta_{lms}\right]$
$\displaystyle+\Delta^{gel}_{+}N^{gel}_{+}+\Delta^{gel}_{-}N^{gel}_{-}+\Delta^{gel}_{0}N^{gel}_{0}$
(1)

where $\beta=1/k_{B}T$ is inverse thermal energy, $\phi_{lms}$ is the volume
fraction of rank $lms$ clusters, $N_{lms}$ is the number of rank $lms$
clusters, $\Delta_{lms}$ is the free energy of formation for clusters of rank
$lms$, and $\Delta^{gel}_{i}$ is the free energy change of species $i$ upon
association to the gel. The free energy in Eq. (1) contains just three
essential pieces of physics: the ideal entropy of mixing for the distribution
of clusters (first term), the free energy of formation of the clusters (second
term), and the free required for species to join the gel (last three terms).

The WiSE is modelled to reside on a lattice containing $\Omega$ lattice sites,
with volume, $v_{0}$ (equivalent to the molecular volume of a water molecule).
Anions and cations each require $\xi_{+}=v_{+}/v_{0}$ and
$\xi_{-}=v_{-}/v_{0}$ lattice sites per cation or anion, respectively, where
$v_{i}$ is the molecular volume of species $i$. Thus, the volume fraction of a
rank $lms$ cluster can be written simply as
$\phi_{lms}=(\xi_{+}l+\xi_{-}m+s)N_{lms}/\Omega$.

The ions and water molecules are modelled to form associations via a fixed
number of association sites, known as the species functionality, $f_{i}$. We
model the cations to have a functionality of $f_{+}$, and anions to have a
functionality of $f_{-}$. We assume that water molecules coordinate to cations
with a functionality of 1 ($f_{0}=1$). The clusters that are formed are
assumed to be Cayley trees (branched with no intra-cluster loops), which are
decorated with water molecules residing on the open cationic association
sites. Thus, the backbone of the cluster purely comprises of ions. Any cluster
containing $l$ cations and $m$ anions will contain exactly $l+m-1$ ionic
associations. As we assume the cations here to be “sticky”, the number of
water molecules in a cluster, $s$, is thus the total number of cationic
association sites in the cluster ($f_{+}l$) minus the number of ionic
associations ($l+m-1$), i.e., $s=f_{+}l-l-m+1$. Note that because
$s=f_{+}l-l-m+1$, the rank of the cluster, $lms$, can actually be specified by
just the ion indices, $lm$. However, for book-keeping reasons, we will
continue to use $lms$ to specify the cluster index.

The equilibrium distribution of clusters can be derived by enforcing chemical
equilibria (shown explicitly in the SI), yielding the following relation:

$\displaystyle\phi_{lms}=K_{lms}\phi_{100}^{l}\phi_{010}^{m}\phi_{001}^{s}$
(2)

where $K_{lms}=\exp(l+m+s-1-\Delta_{lms})$ is the cluster equilibrium
constant. Thus, the equilibrium cluster distribution depends on the explicit
form of $\Delta_{lms}$, which is derived in the SI. The strength of the
presented theory comes from the fact that $\Delta_{lms}$ only requires six,
physically transparent parameters [$\xi_{\pm}$, $f_{\pm}$, $\tilde{\lambda}$
(the ionic association constant) and $\tilde{c}_{salt}$ (dimensionless
concentration of salt - $\\#$ per lattice site)] to be specified to determine
all possible $\Delta_{lms}$.

Provided the functionalities, $f_{\pm}$, of anions and cations are greater
than two, the clusters can become infinitely large. For this to occur, the
probability that a cation (anion) has another cation (anion) separated by an
anion (cation) must reach 1. The point at which this occurs (i.e. the gelation
point) is given by the condition

$\displaystyle p^{*}_{+-}p^{*}_{-+}(f_{+}-1)(f_{-}-1)=1$ (3)

where $p^{*}_{+-}$ and $p^{*}_{-+}$ are the critical (denoted with $*$) ion
association probabilities for gelation. The ion association probability,
$p_{ij}$ is defined as the fraction of association sites on species $i$
occupied by an association to species $j$. Equation (3) is a classical result
of Flory Flory (1941) and Stockmayer Stockmayer (1943) and is analogous to
bond percolation on a Bethe lattice which arises when the cumulative
probability of a cluster to branch indefinitely becomes unity. These
association probabilities essentially close the theory, and make the
connection between microscopic bonding picture and bulk thermodynamics. As
shown in the SI, these association probabilities can be explicitly expressed
as a function of salt concentration:

$\displaystyle
f_{+}p_{+-}=f_{-}p_{-+}=\frac{1+\mathcal{D}\tilde{c}_{salt}-\sqrt{(1+\mathcal{D}\tilde{c}_{salt})^{2}+4f_{-}\tilde{\lambda}\tilde{c}_{salt}}}{2\tilde{c}_{salt}(\tilde{\lambda}-1)}$
(4)

where
$\mathcal{D}=f_{+}(\tilde{\lambda}-1)+f_{-}\tilde{\lambda}-\xi_{+}-\xi_{-}$
and the ionic association constant

$\displaystyle\tilde{\lambda}=\exp\left\\{-\beta(\Delta U_{+-}-\Delta
U_{+0})\right\\}$ (5)

governs the extent of ion association, where $\Delta U_{+-}$ is the energy of
a cation-anion association and $\Delta U_{+0}$ is the energy of a cation-water
association. The ion association constant is determined by the difference
between of $\Delta U_{+-}$ and $\Delta U_{+0}$ as a direct result of the
“sticky” cation approximation: every time a cation-anion association is
formed, a cation-water association must be destroyed. The cation-anion and
cation-water associations are largely electrostatic owing to the monopolar and
dipolar nature of the ions and water molecules, respectively.

Also as a consequence of the “sticky” cation approximation, the probabilities
corresponding to cation-water associations, $p_{+0}$ and $p_{0+}$, can be
directly related to the ion association probabilities:

$\displaystyle
p_{+0}=\left(\frac{\tilde{c}_{water}}{f_{+}\tilde{c}_{salt}}\right)p_{0+}=1-p_{+-}$
(6)

where the dimensionless concentration of water is given by
$\tilde{c}_{water}=1-(\xi_{+}+\xi_{-})\tilde{c}_{salt}$, as prescribed by our
assumption of electrolyte incompressibility.

In WiSEs, the salt concentration is readily altered through the
addition/removal of salt or water. At low salt concentrations, one does not
expect to form an ionic gel, but in highly concentrated WiSEs that are
employed in battery applications, a gel could potentially form. From solving
Eq. (3), we can obtain the critical salt concentration,
$\tilde{c}^{*}_{salt}$, required to form an ionic gel:

$\displaystyle\tilde{c}^{*}_{salt}=(\mathcal{A}\tilde{\lambda}+\mathcal{B})^{-1}$
(7)

where $\mathcal{A}$ and $\mathcal{B}$ are constants given by

$\displaystyle\mathcal{A}=\frac{1}{2}\sqrt{\frac{f_{+}f_{-}}{(f_{+}-1)(f_{-}-1)}}+2\sqrt{f_{+}f_{-}(f_{+}-1)(f_{-}-1)}-f_{+}-f_{-}$
(8)

and

$\displaystyle\mathcal{B}=\xi_{+}+\xi_{-}-\frac{1}{2}\sqrt{\frac{f_{+}f_{-}}{(f_{+}-1)(f_{-}-1)}}$
(9)

Alternatively, one can think about fixing the concentration of salt in the
WiSE and varying the temperature. One can invert Eq. 7 and obtain a critical
value of $\tilde{\lambda}^{*}=-\mathcal{B}+1/\mathcal{A}\tilde{c}_{salt}$,
which would directly correspond to a critical gel-point temperature for a WiSE
at a specific salt concentration:

$\displaystyle T*=\left(\frac{\Delta U_{+-}-\Delta
U_{+0}}{k_{B}}\right)(\mathcal{B}+1/\mathcal{A}\tilde{c}_{salt})^{-1}$ (10)

This equation can be used to predict the critical gel line in the
concentration–temperature plane. However, we must be sure to take extreme care
to ensure that the sticky cation approximation is applicable across the
studied temperature range. At high temperatures, it becomes more likely that
association sites on the cation are unoccupied (with either anion or water),
which could have a considerable effect on the predicted critical gel line.

When the gel is present, it becomes necessary to determine the exact fraction
of species in the gel, as well as the distribution of finite clusters not
incorporated in the gel (i.e. sol). This computation requires an slightly
altered procedure than what was outlined in this section. For the details
outlining this procedure, as well as an extended discussion of the critical
gel line in the concentration–temperature plane, we direct the reader to the
SI.

##  III Results

The theory presented in the preceding section allows us to predict the ion
clustering and solvation in WiSEs. In order to test the validity of our
developed theory, we performed a series of molecular dynamics (MD) simulations
of aqueous lithium bis(trifluoromethanesulfonyl)imide (LiTFSI), lithium
triflouromethanesulfonate (LiOTF), and lithium bis(fluorosulfonyl)imide
(LiFSI) electrolytes for concentrations ranging from the dilute to “water-in-
salt” regime. In this section, we summarize our results and observations from
these MD simulations, and then we compare those results against our theory. We
mainly focus on the results of LiTFSI in the main text, to avoid repeating
similar observations in the other salts. The results for LiOTF and LiFSI
solutions are mainly shown in the SI.

###  III.1 Criteria and Characteristics of Association

In order to study the degree of ionic aggregation and solvation in these
electrolytes we first need a precise definition of what constitutes a lithium-
water association and a lithium-anion association. The Li+ interacts very
strongly with the oxygen of the water molecule, as well as the oxygens of the
TFSI-, OTF-, and FSI- ions. The affinity of Li+ to water or anionic oxygens
can be seen clearly by the spatial distribution functions visualized in Fig.
1a as iso-density surfaces. These iso-density surfaces encapsulate regions
around the water molecule or anion where Li$+$ density is on average _at
least_ 4$\times$ its bulk density. It can be seen clearly that these regions
are typically around the oxygens of the TFSI-, OTF-, and FSI- ions. Thus, it
is natural to define a lithium-anion or lithium-water association based on the
proximity of Li+ to oxygen atoms, belonging either to anions or water,
respectively. In particular, we use a Li+-O distance threshold of 2.7 $\AA$ to
determine if a Li+ is associated to oxygens on water molecules or anions.

This distance threshold is not arbitrary, and can be understood clearly when
examining the nearest neighbor distance distributions between the associating
species in the electrolyte. Here we will define an A–B nearest neighbor
distance distribution to be the probability of finding a species A, whose
nearest neighbor, B, is a distance, r, away from A. In Fig. 1b & c, we plot
the Li+–O(TFSI-) (Fig. 1b) and the O(water)–Li+ (Fig. 1c) nearest neighbor
distance distributions for concentrations of LiTFSI ranging from 0.28m–21m.

Figure 1: (a) The spatial distribution functions of Li+ around water, TFSI-,
OTF-, and FSI- showing iso-density surfaces corresponding to 4x the bulk
density of Li+ . The probability distribution of the nearest-neighbor
distances for (b) Li+–O(TFSI-) and (c) O(water)–Li+. The insets in (b) and (c)
zoom in on the probability distributions within the boxed region.

For the Li+–O(TFSI-) nearest neighbor distance distribution first has a sharp
peak at a distance of roughly 2 $\AA$, approximately corresponding to the
distance of closest approach between lithium and oxygen. The distribution then
decreases to a minimum at 2.7 $\AA$, before we observe another peak centered
around distances of roughly 4 $\AA$, which becomes sharper as the
concentration increases. Thus, there exists one sub-population of lithium in
direct contact with anion and another sup-population of lithium not in direct
contact with anion. We have a near identical story for the O(water)–Li+nearest
neighbor distance distribution. The distribution is indeed bi-modal, with two
peaks again corresponding to sub-populations of water molecules that are in
direct contact (2 $\AA$) with lithium and another that is not in direct
contact ($\sim$4 $\AA$) with lithium, separated by probability minimum at
roughly 2.7 $\AA$. Thus, our Li+-O cutoff distance of 2.7 $\AA$ provides a
clear-cut criterion for both Li+-anion and Li+-water associations.

The insets of Fig. 1 b&c zoom in on the behavior of the nearest neighbor
distributions at large distances. We observe that at high concentrations ($>$2
m), the Li+–O(anion) nearest neighbour distance distributions rapidly decay to
zero for nearest-neighbor distances greater than 4$\AA$. Therefore, at these
high concentrations, it becomes exceedingly likely that lithium is present
only in one of two states: in direct contact with an anionic oxygen (2 $\AA$),
and indirect second-shell contact with an anionic oxygen (4 $\AA$). However,
at low salt concentrations ($<$2 m), the Li+–O(anion) nearest-neighbor
distance distributions do not rapidly decay after the second maximum at 4
$\AA$. Instead, they actually observe a third peak in the nearest neighbor
distribution, which is much more spread out. Thus, it is clear that at
concentrations less than 2 m there are 3 characteristic states in which ions
exist: direct first-shell contact (2 $\AA$), indirect second-shell (i.e.
solvent-separated Eigen and Tamm (1962); Marcus and Hefter (2006); Seo  _et
al._ (2012a, b)) contact (4 $\AA$), and free ($>$6 $\AA$). In all cases, the
probability for any lithium ion to have its nearest anionic oxygen be very
large distances ($>12\AA$) away becomes extremely small.

We should note that there is no distinction in our theoretical model between
ions that are in indirect second-shell contact (4 $\AA$) and completely free
ions($>$6 $\AA$); these scenarios would both correspond to ions that are not
considered to be associated. The exchange between second-shell and free states
could likely be modelled well within Bjerrum’s framework of ion pairing
Bjerrum (1926), but it would only be useful to model this exchange at low salt
concentrations. For concentrations greater than 2 m, the probability of any
ions being significantly separated ($>$6 $\AA$) becomes exceedingly rare, and
the majority of ions exist only in direct first-shell contact or indirect
second-shell contact. For these high salt concentrations, our binary
definition of ion association becomes much more appropriate; nearest ionic
neighbors are either in the first shell or the second shell.

In Fig. 1e, there is a similar story for the O(water)–Li+ nearest-neighbor
distance distributions. At all concentrations there are peaks corresponding to
water molecules in the first (2 $\AA$) and second solvation (4 $\AA$) shells
of Li+. However, for concentrations less than 2 m, (as seen in the inset of
Fig. 1 c), we observe a peak in the O(water)–Li+ nearest neighbor distribution
corresponding to water molecules not within the the first two coordination
shells of any Li+, i.e. truly free water. Again our model treats water
molecules outside of the first coordination shell as not being associated to
Li+ and does not distinguish between second shell water molecules and truly
free water molecules (outside of the second solvation shell of all Li+). There
are existing statistical models based on multi-layer Langmuir adsorption
Dutcher  _et al._ (2011) that can handle multi-shell ion hydration. However,
such an extension to our model would only be necessary for low salt
concentrations, as it becomes exceedingly unlikely for water molecules to
exist outside the second shell of Li+ for concentrations greater than 2 m.

Thus, our simulations show that our model’s simple binary picture of species
as associated or not associated (either to ion or water) is not strictly true
for low salt concentrations. However, this binary picture does become valid at
high salt concentrations, in which species reside almost strictly two distinct
states corresponding to species that are associated or unassociated. The
primary focus of this paper is on super-concentrated WiSE’s, where our model
is likely to be largely valid, as indicated by our simulations.

###  III.2 Li+ Coordination: Ionic Association vs. Solvation

A very important phenomenon in WiSEs is the competition between ionic
association and solvation. From our simulations, we find that Li+ ions
generally have four associations to oxygen atoms; either to the anions or
water molecules, resulting in a competition between Li+–O(anion) and Li+–OH2
coordination due to the limited space around the Li+ ions. If Li+–OH2
associations dominate over Li+–O(TFSI-), then ion aggregation is suppressed
and the formation of an ionic gel is forestalled. The competitive coordination
can be seen clearly in Fig. 2a, where we have measured (via MD simulation) the
average coordination number of Li+ to the water oxygen or TFSI- oxygens as a
function of the salt concentration. We observe that the average Li+–OH2 and
Li+–TFSI- coordination numbers respectively decrease and increase
monotonically as a function of salt concentration. This is further illustrated
in Fig. 2b & c, where the Li+ coordination number probability distributions
are plotted for both water oxygens (Fig. 2b) and TFSI- oxygens (Fig. 2c). For
any given concentration, the probability distributions for Li+ coordination by
anion and water have a strong reflective symmetry with each other, which is a
result of the approximately conserved Li+-O coordination number of four. Note,
that this Li+-O coordination number conservation is in strong agreement with
the “sticky” cation approximation that was a main assumption of the derived
theory.

Figure 2: (a) The average coordination number of Li+by water oxygens and by
TFSI- oxygens as a function of concentration, showing the competitive nature
of ion association and solvation. (b) The probability distributions for the
number of TFSI- oxygens coordinating lithium. (c) The probability
distributions for the number of water oxygens coordinating lithium.
Distributions are plotted for LiTFSI concentrations of 0.28 m-21 m. (d) A
schematic showing cartoon depictions of typical coordination environments of
Li+ in which Li+ is coordinated by 4 oxygen atoms either from TFSI- or water.

This concept is further depicted in Fig. 2d, in which we draw the 5 possible
coordination environments that conserve 4 associations to Li+. Note that Li+
coordination is not the complete picture. Each TFSI- ion coordinating Li+
contains additional oxygen atoms that can associate with further Li+. On the
other hand, each of the hydrating water molecules serve as a network
terminating agent because they can only associate to a single Li+ ion. Thus,
Li+ coordination informs us about the average number and types of associations
that Li+ makes, but it does not tell us the extent and size of clusters that
are formed.

###  III.3 Ionic Aggregation and Gelation

As was previously mentioned, the fact that Li+ and TFSI- (as well as OTF- and
FSI-), are able to form multiple associations permits the formation of large,
potentially infinite, ion networks. In Fig. 3a we plot the probability
distribution, $\mathcal{P}_{n}$, for an ion to be incorporated in an ion
cluster containing a total of $n$ ions (anions + cations), for a range of
different LiTFSI concentrations. We see that for concentrations below 19 m,
$\mathcal{P}_{n}$ decreases monotonically as a function of the cluster size.
This observation might be intuitively expected as larger clusters require the
coordination of a larger number of ions, which would seem increasingly
unlikely at low concentrations. For increasing salt concentrations (less than
19 m), it becomes increasingly more likely to observe higher order clusters.
This is intuitive because ions are more likely to associate when there is more
salt and less water in the mixture.

At concentrations above 19 m, the cluster size distributions are no longer
monotonically decreasing. Rather, we observe a remarkable peak in the
distribution corresponding to extremely large ionic clusters. This observation
is illustrated more clearly by the cumulative distribution function (CDF)
plotted in Fig. 3b, where we see a seemingly anomalous uptick in the CDF
indicated by the dashed line enclosing the shaded region. In fact, the
clusters within the probability distribution peaks (indicated by the gel
label), as well as the shaded region of the CDF, correspond to clusters that
involve a large fraction of the total amount of ions in the simulation. We can
see this in Fig. 3c, where the largest cluster is highlighted in blue in a
snapshot of a MD simulation of 21 m LiTFSI. The highlighted cluster extends
throughout the entire simulation box, and we would expect it to grow
arbitrarily large depending on the size of the simulation. Thus, we define
this percolating structure as the _ionic gel_. This graphic (in Fig. 3c) also
illustrates that the ionic gel is far from crystalline. In fact, the gel is
quite spatially disordered, which is characteristic of polymeric gels Flory
(1953). This spatial disorder comes through the chaotopic anions, not the
kosmotropic cations. Physically, what is propelling the percolation of the ion
cluster is that ion associations become so abundant in the mixture that
clusters have a finite probability to continue indefinitely as opposed to
terminating.

Figure 3: (a) The probability distribution of ionic cluster sizes for a range
of different concentrations (0.2 m-21 m) of LiTFSI. Peaks in the distributions
for concentrations greater than 16 m correspond to percolating ionic gels. (b)
The cumulative distribution function (CDF) of ionic cluster sizes for the
various concentrations. The shaded region corresponds to ionic gel which
result in a sharp uptick in the CDF. (c) A visualization of an MD snapshot of
21 m LiTFSI. The ionic gel is highlighted in blue, and is observed to
percolate throughout the entire simulation box.

In polymeric mixtures, the presence of gel can lead to large changes in the
mixture’s physical properties, such as the divergence of viscosity and the
onset of finite elasticity Flory (1953). In WiSE’s the ionic gels are patched
together by attractive electrostatic interactions, which are thermally
reversible (breakable). Thus, we hypothesize that WiSE’s containing gels are
potentially viscoelastic fluids, perhaps able to bridge the gap between
typical liquid electrolytes and polymeric solid-state electrolytes in a
tunable fashion. However, the viscoelastic properties of WiSEs would depend
highly on the lifetime and strength of the associations. For example, Winter
Winter (2002) describe such physical gels as having solid-like behavior below
a certain yield stress (related to the strength of associations) and at time
scales shorter than the so-called renewal time of the gel (related to the mean
association lifetime). We find that the association lifetimes are on the order
of picoseconds (3-13ps for aqueous LiTFSI electrolyte as shown in figure S4 in
the SI). The reason for the relatively short association lifetimes is not
because the associations are weak, but rather because there is a strong
competition between ion association and lithium hydration, as we have
explained, which tends to limit the lifetime of the ion associations.
Actually, as shown in the subsequent section, the success of the sticky-cation
approximation in describing the clustering in lithium-based WiSEs implies that
ion associations are very strong, with magnitudes well above $k_{B}T$.
Nonetheless, the picosecond scale lifetimes suggest that ion network-forming
WiSEs could display solid like properties, but only below picosecond
timescales.

Perhaps of equal importance is what the association lifetime tells us about
the ability of ion clusters to contribute to diffusion in our system. In ref.
44, it was shown that species contribute to diffusion only if their lifetimes
are longer than the decay of their velocity auto-correlation functions
(VACFs). In ref. 42, this idea was extended to ion clusters yielding an
approximate constraint on the sizes of clusters that may contribute to the
diffusive relaxation of the system:

$(l+m)^{2/3}(l+m-1)<\tau_{B}/\tau_{\nu}$ (11)

where $l$ is the number of cations in the cluster, $m$ is the number of anions
in the cluster, $\tau_{B}$ is the association lifetime, and $\tau_{\nu}$ is
the characteristic time of the decay of the VACF for free ions (generally
around 1 ps). Thus, for ion association times of 3-13 ps (as is the case for
aqueous LiTFSI electrolyte), we would expect only clusters containing 3-5 ions
_or less_ could contribute diffusively to the system.

###  III.4 Comparison of Theory and Simulation

Our MD simulations demonstrate the relevance of ionic gelation, and
competitive ionic association and solvation in WiSEs. Our thermodynamic model
was designed to describe such phenomena. Thus, in this section, we test the
ability of our model to capture our observations from the aforementioned
simulations.

Our model contains a number of molecular parameters: $\xi_{\pm}$, $f_{\pm}$,
$\tilde{\lambda}$. The size parameters $\xi_{\pm}$ are obtained by computing
the molecular volumes consistent with the force fields in our MD simulations
(volume enveloped by the overlapping Lennard-Jones radii of the molecule). The
Li+ functionality was determined to be 4, which is consistent with the Li+
coordination data in Fig. 2. Although we did observe rare instances of Li+
coordinating more than 4 species (water+anionic oxygens), these instances were
rare enough to still model Li+ as having $f_{+}=4$. The anionic functionality
can be reasoned clearly by the observed SDFs in Fig. 1a. Here we can see the
regions around each of the anions with a specific preference for Li+. For each
anion, (TFSI-, OTF-, and FSI-), we observed 3 distinct regions, indicating
that they will have all have an ion functionality of $f_{-}=3$. This is
slightly unexpected for TFSI- and FSI-, because these anions have 4 oxygens
that might associate to Li+, which would ostensibly indicate that $f_{-}=4$
for these species. Nevertheless, it seems to be the case that only 3 of the
anionic oxygens can bind an Li+ at any one time. Again, there are rare
instances where the anions associate to more than 3 Li+ ions, but these
instances are exceedingly rare. The functionality of the solvent, $f_{0}$ was
determined to be 1.

Figure 4: Association probabilities ($p_{+-}$, $p_{-+}$, $p_{+0}$, & $p_{0+}$)
are plotted as a function of salt concentration for (a) LiTFSI (b) LiOTF and
(c) LiFSI electrolytes. The square data points are obtained from molecular
dynamics (MD) simulations, and the solid curves are generated from our
theoretical model, using association constants, $\tilde{\lambda}$, that are
computed from the MD simulations. All curves were generated using the
parameters, $\xi_{+}=0.4$, $f_{+}=4$, and $f_{-}=3$. LiTFSI curves used the
additional parameters $\lambda=0.075$ and $\xi_{-}=10.8$. LiOTF curves used
the additional parameters $\lambda=0.184$ and $\xi_{-}=5.9$. LiFSI curves used
the additional parameters $\lambda=0.049$ and $\xi_{-}=6.8$.

The remaining parameter, $\tilde{\lambda}$ can also be easily computed from MD
simulation data. As is shown in the SI, the ion association constant is
related to the association probabilities, $p_{ij}$, through a mass action law
which describes the exchange of water and anions in the lithium coordination
shell:

$\displaystyle\tilde{\lambda}=\frac{p_{-+}(1-p_{0+})}{p_{0+}(1-p_{+-})}$ (12)

Thus, in order to compute $\tilde{\lambda}$ we need only compute the
association probabilities from simulation. The association probabilities can
be computes simply using the following formulae:

$\displaystyle p_{ij}=\left\langle\frac{\\#\text{ of associations of type
}ij}{f_{i}\cdot\\#\text{ of molecules of type }i}\right\rangle$ (13)

where the $\langle\cdot\rangle$ operation indicates an ensemble average. The
probabilities are obtained as a function of salt concentration, and then
$\tilde{\lambda}$ is computed from Eq. (12) and averaged across all
concentrations.

In Fig. 4, we plot the association probabilities, $p_{ij}$, as a function of
salt concentration for each of the simulated electrolytes, LiTFSI, LiOTF, and
LiFSI. For each salt, the theory curves (with model parameters computed from
MD) match the simulated values almost quantitatively in the entire
concentration range studied.

With the parameters in our model fully specified, we may explore theoretical
predictions of the model and test their accuracy with MD. Of foremost
importance, is the clustering of ions as a function of salt concentration. In
Fig. 5a, we plot the concentration dependence of the fractions of various ion
cluster types: free ions (not bound to any other ions), ion pairs (one cation
bound to one anion), finite clusters (clusters larger than ion pairs up to,
but excluding, the gel) and ionic gel (as previously defined). The model does
a remarkable job reproducing the trends we observe from simulations for all of
the plotted ion cluster types.

As expected for aqueous electrolytes in dilute to moderate salt concentrations
($<$5 m), the electrolyte is dominated by free ions, and the ion aggregates
that appear are almost entirely ion pairs. In this regime, ion aggregation
would be consistent with Bjerrum’s description of ion pairing. However, upon
increasing the salt concentration to $\sim$10 m, we observe that ions are more
likely to appear in high order clusters than simple ion pairs. Though, the
free ions still comprise the majority of the total ions in the electrolyte.
Further, increasing the concentration up to roughly 19 m, we find that ions
are more likely to be present in high order finite clusters than both ion
pairs and free ions. In fact, the fraction of ions appearing in ion pairs
begins to decay, with the max fraction of ion pairs occurring at a
concentration of roughly 12 m. Finally, for concentrations greater than
roughly 19 m, we enter the ionic gel regime. The gel point is characterized in
Fig. 5a not only by the onset of a finite gel fraction, but also by a sharp
decrease in the fraction of ionic cluster larger than ion pairs, suggesting
that large portions of the high order finite aggregates combine to form the
gel at and beyond the gel point.

The weight-averaged degree of ionic aggregation, $\bar{n}_{w}$, which is
plotted as a function of salt concentration in Fig. 5b, also provides insight
into the extent of clustering. The weight average degree of aggregation is
defined as the following

$\displaystyle\bar{n}_{w}=\langle
l+m\rangle=\frac{\sum_{lm}(l+m)\alpha_{lm}}{\sum_{lms}\alpha_{lm}}$ (14)

where $\alpha_{lm}$ is the fraction (probability) of ions incorporated in
ionic clusters containing $l$ cations and $m$ anions (see SI for details). The
phrase “weight-averaged” arises because $\alpha_{lm}$ is weighted by the total
number of ions in the cluster. Of course $\alpha_{lm}$ and thus $\bar{n}_{w}$
may be computed from MD directly. The theoretical weight-averaged degree of
aggregation is derived in the SI, and can be written analytically in closed-
form, first derived by Stockmayer in 1952 Stockmayer (1952). Note that in Fig.
5b, we are plotting the weight-averaged degree of average of only the _finite_
aggregates in the sol. Thus for the MD measurements, we exclude the largest
cluster at each timestep when computing $\alpha_{lm}$. We observe that the
simulated weight-average degree of aggregation has a maximum occurring at 19
m, which is in good agreement with the theoretically obtained value of 19.5 m.
Moreover, the theoretical and simulated weight-average degrees of aggregation
agree very well with each other across the entire concentration range. Though,
the MD measurements do not display true divergence at the gel point, like the
theoretical prediction, due to the finite size of our simulation Christensen
and Moloney (2005). This finite simulation box effectively caps the size that
finite ion clusters, preventing the weight-averaged degree of aggregation from
truly diverging.

Figure 5: (a) The fractions of various ionic cluster types (free ions, ion
pairs, finite high order clusters and ionic gel), (b) the weight-average
degree of aggregation (for finite clusters), and (c) the coordination numbers
of Li+ by oxygen atoms from water and oxygen atoms by TFSI-,. The curves were
generated using the parameters $\lambda=0.075$, $\xi_{+}=0.4$, $\xi_{-}=11.3$,
$f_{+}=4$, and $f_{-}=3$.

In Fig. 5c we plot the Li$+$ coordination by water (blue curve) and the Li$+$
coordination by anion (red curve) as functions of salt concentration, and we
observe a quantitative match of the theory and MD simulations. As we
previously mentioned, these curves demonstrate the competition between ion
solvation and ion association. The inclusion of this essential physics in our
model is responsible for a large majority of its success. It also should be
noted that the onset of ionic gelation does not leave any characteristic
signature in the coordination numbers. Instead, we observe the coordination
numbers of Li$+$ by both water and anion progress smoothly through the
critical gel point concentration.

##  IV Discussion

###  IV.1 Origins of Charged Nano-Domains in WiSEs

The existence of percolating ionic networks in WiSEs has been noted in
literature. For example, in Refs. 26; 71; 28; 29, the authors performed MD
simulations and found percolating ion networks in super-concentrated
electrolytes when the ions contained so-called khaotropic anions, such as
TFSI-. Moreover, in Ref. 18, vibrational spectroscopy indicated the existence
of high order aggregates and inter-connected ion networks in super-
concentrated LiTFSI electrolytes. It is clear that our model and simulations
are congruent with these observations of percolating ion networks via its
prediction of gelation in WiSEs at high enough salt concentrations.

However, in addition to the percolating ion networks, Refs. 19; 27; 72; 73 all
discuss the emergence of lithium and water-rich nano-domains, which allow for
the facile transport of Li+ within the electrolyte. Interestingly, our model
can shed light on the formation of these nano-domains, as well. Within our
theoretical framework, the electrolyte species that are not associated to the
gel are referred to as the _sol_ , which consists of just free species and
finite clusters. We can think about the sol as potentially forming a second
network in the electrolyte, if the sol spatially percolates through the
electrolyte. In this case, the sol would form a network not via ionic
associations, but rather by the virtue of simply occupying enough volume to
ensure spatial percolation through the electrolyte.

To determine if the sol network spatially percolates, a “site” percolation
model can be used, where the electrolyte is divided into lattice-like voxels
with some probability of the voxel being occupied by the sol. This probability
of a voxel being occupied by the sol is equivalent to the volume fraction of
the sol, $\phi_{sol}$ ($\phi_{sol}=1-\phi_{gel}$). It is clear that when
$\phi_{sol}$ is large enough (or conversely $\phi_{gel}$ is small enough), the
sol will easily percolate through the electrolyte. For example, for salt
concentrations below the gel point, the entire electrolyte is sol, and it
trivially percolates through the electrolyte. However, when $\phi_{sol}$ falls
below a certain threshold, site percolation theory predicts that te sol would
no longer percolate through the electrolyte. The critical $\phi_{sol}$ for
site percolation is not certain, as our electrolyte does not reside on a
crystalline lattice structure. However, critical values for simple lattice
structures will range from $\phi^{*}_{sol}=0.2-0.31$ (Simple cubic: 0.31, BCC:
0.25, and FCC: 0.20) Christensen and Moloney (2005). We note that the
percolation of sol could also be framed in the context of continuum
percolation theory, which would not require the specification of a lattice,
but rather the shapes of the percolating species. Still, continuum percolation
in three dimensions would yield roughly the same critical percolation
thresholds (between 0.19 and 0.3 for spheroids with aspect ratio less than
four Garboczi  _et al._ (1995)).

Interestingly, this sol “network” will share many of the same properties of
the water and lithium rich nano-domains identified in Refs. 19; 27; 72; 28;
73.First, when $f_{+}>f_{-}$, the sol will have a net positive charge (as
noted in Refs. 19; 28) containing more Li+ than anion. Second, the sol will
also contain increasingly _free_ ionic species as the ionic gel becomes more
substantial McEldrew  _et al._ (2020b). Lastly, when these free species are
majority Li+ (as will be the case for $f_{+}>f_{-}$), the model ensures that
the free Li+ is fully hydrated by water. These properties, yield a sol network
that is water and Li+ rich when the ion gel network is substantial and
$f_{+}>f_{-}$. However, in order for this water and Li+ rich network to
facilitate fast transport of Li+ it must _percolate_ through the electrolyte.
Otherwise, the sol will exist only in disconnected pockets, and ion transport
would likely require more activated processes such as association/dissociation
from the gel.

In Fig. 6, we show several important concentration regimes emerging in LiTFSI.
The most dilute region is the Salt-in-Water Electrolyte (SiWE) regime, where
the volume fraction of the salt is less than 0.5 (the majority of the volume
of the electrolyte is water). For concentrations ranging from 5-20 m, we have
a WiSE without a percolating ion gel (the sol will trivially percolate through
the electrolyte in this regime). Upon increasing the salt concentration past
the critical gel concentration ($\sim$20 m), an ionic gel emerges and grows
rapidly in mass and volume. From about 20 m-27.5 m the WiSE exists with both
percolating ion networks, as well as percolating sol. For concentrations
greater than about 27.5 m, the volume fraction of sol ($1-\phi_{gel}$) reduces
below 0.25 (BCC lattice assumed here) and the sol no longer percolates through
the electrolyte. In this regime the transport properties of the WiSE could
substantially altered, as the pockets of sol will become trapped within the
strongly gelled electrolyte. For example, when the sol is able to spatially
percolate, we can expect that conduction largely occurs via sol-mediated
vehicular transport of ions. When the sol no longer spatially percolates
through the electrolyte, ion conduction would likely depend on alternative
mechanisms, such as activated hopping mechanismsAndersson  _et al._ (2020) or
ion-exchange between sol and gel phases, similar to the “ionic semiconductor”
concept introduced in ref. 44.

Figure 6: The volume fraction of gel (blue curve), sol (red curve) and total
salt (yellow curve) in LiTFSI electrolyte are plotted as a function of salt
concentration. Various concentration regimes are depicted. The Salt-in-Water
Electrolyte (SiWE) regime spans volume fractions of salt less than 0.5. The
Water-in-Salt Electrolyte (WiSE) pre-gel regime spans salt volume fractions
greater than 0.5 but concentrations below the critical gel concentration. The
WiSE gel regime spans WiSEs with concentrations above the critical gel
concentration. WiSEs with a sol volume fraction greater than 0.25 (for BCC
lattice) will have a percolating sol phase. WiSEs with a sol volume fraction
of less than 0.25 will not have a percolating sol phase. The theoretical
curves were generated using the parameters $\lambda=0.075$, $\xi_{+}=0.4$,
$\xi_{-}=11.3$, $f_{+}=4$, and $f_{-}=3$.

###  IV.2 Ion transference number at super-concentration

It is well known that the ionic conductivity of WiSEs is strongly non-
monotonic, peaking at intermediate concentrations, before strongly decreasing
at higher concentrations Suo  _et al._ (2015); Ding  _et al._ (2017); Borodin
_et al._ (2017a); Ding and Xu (2018); Horwitz  _et al._ (2020). The reason for
this non-monotonicity has been attributed to the increasing viscosity at high
salt concentrations, as well as the substantial ionic clustering in WiSEs,
which reduces the “mobile” charge carriers in solution. However, what is less
well established is the asymmetry of the ionic transference numbers in WiSEs.
It has been observed in LiTFSI electrolytes that the Li+ transference number
increases as a function of concentration. This asymmetry has been attributed
to the emergence of asymmetric ionic clusters, resulting in an unequal
dissociation of Li+ and TFSI-, favoring more free Li+ Borodin  _et al._
(2017a); Yu  _et al._ (2020); Horwitz  _et al._ (2020).

As eluded to in the previous section, our model is able to rationalize the
presence of asymmetric clusters in WiSEs, which had not been established
clearly in literature. If the cations and anions are able to make a different
numbers of associations ($f_{+}\neq f_{-}$), then the associations will be
asymmetrically distributed amongst them. Specifically, when $f_{+}>f_{-}$,
cations will necessarily have a lower ion association probability than anions
($p_{+-}<p_{-+}$), which will inevitable yield a greater number of free
cations than anions. This asymmetry is exacerbated when the ionic gel is
present, as we described in the previous section. The unequal ion
functionalities result in the gel and sol phases with a non-zero, but equal
and opposite, net charge. For $f_{+}>f_{-}$, this results in a net positive
sol, and a net negative gel, as has been observed in literature Borodin  _et
al._ (2017a); Lim  _et al._ (2018); Zheng  _et al._ (2018); Yu  _et al._
(2020); Gonzales  _et al._ (2020); Andersson  _et al._ (2020). Moreover,
because the ionic gel is percolating and macromolecular, we can assume that it
does not have translational degrees of freedom and thus will not contribute to
the ionic current. On the other hand, the net positive sol will dominate the
ionic current, and yield cationic transference numbers greater than 0.5
McEldrew  _et al._ (2020b, a). Furthermore, the ion association lifetimes
(3-13 ps) indicated that only clusters containing less than 3-5 ions (as
estimated by eq. (11)) live long enough to conduct current, with free ions
dominating the population of those current-conducting species. When free ions
dominate the conductivity Feng  _et al._ (2019); McEldrew  _et al._ (2020b,
a), we can approximate the ionic transference number in terms of the free
cation and anion fractions, $\alpha_{\pm}$, and the cation/anion self-
diffusion coefficient ratio, $\mathcal{R}=D_{+}^{0}/D_{-}^{0}$:

$\displaystyle
t_{+}=\frac{\alpha_{+}\mathcal{R}}{\alpha_{+}\mathcal{R}+\alpha_{-}}.$ (15)

This formula can be used to estimate the transference number in WiSEs using
computed free ion fractions from theory or from MD simulations, and by knowing
the self-diffusion coefficient ratio in the dilute regime, where the free ion
fractions of anions and cations will be approximately equivalent.

Of course the cation transference number may be computed explicitly from the
MD simulation via the corrresponding Green-Kubo formula or equivalent mean-
square displacement. However, these computations have been performed elsewhere
Borodin  _et al._ (2017a); Gonzales  _et al._ (2020). Thus, for simplicity we
use Eq. (15) with free ion fractions as computed from the theory or MD
simulations (note that some more sophisticated equations were derived in Ref.
42). The ratio, $\mathcal{R}$ was specified from the experimental data to be
$1.17$ by examining the experimental transference number at the lowest salt
concentration (where the free cation and anion fractions will be roughly
equivalent), as seen in Fig. 7. Note, the experimental data is obtained from
pulsed field gradient NMR, which does not account for cross-correlations
between ions, thus the transference number is denoted as the “apparent”
transference number. Again the MD simulation curve in fig. 7 uses eq. (15),
but with free ion fractions computed from MD simulation. Nonetheless, we can
see in Fig. 7 that the MD estimated $t_{+}$ is very close to the experimental
values, which strongly implies that the current conducted predominantly by
free species in WiSEs. The theoretical curve qualitatively agrees with the MD
and experimental data, but there is appreciable difference between the
theoretically predicted transference numbers. This is because the theory
actually underpredicts the asymmetry of free cations to free anions.
Nonetheless, the qualitative match of the theory suggests uneven ion
functionalities are extremely important factors in manipulating the ion
transference number.

Figure 7: The “apparent” cationic transference number is plotted as a function
of salt concentration for aqueous LiTFSI electrolyte. The theoretical curve
(smooth line) and MD data (squares connected by dashed line) were generated
from Eq. (15), using $\mathcal{R}=1.17$, as determined from the experimental
data (circles connected by dotted line) from Ref. 77. The theoretical curve
was generated using the parameters $\lambda=0.075$, $\xi_{+}=0.4$,
$\xi_{-}=11.3$, $f_{+}=4$, and $f_{-}=3$.

###  IV.3 Implications for the Electrochemical Stability Window of WiSEs

Perhaps the most exciting property of WiSEs, and super-concentrated
electrolytes in general, is there expanded electrochemical stability window
(ESW). Two primary aspects that govern the ESW of WiSEs is the ability of the
WiSE to form a solid-electrolyte interface (SEI) at the negative electrode,
and the reduced thermodynamic activity of water in the electrolyte. The
formation of the SEI kinetically suppresses the hydrogen evolution at the
negative electrode, and the reduced activity of water decreases the
thermodynamic driving force for the oxygen evolution on the positive
electrode.

The most important aspect of ESWs in WiSEs is its ability to form an SEI
layer. In WiSEs, the electrolyte can no longer rely on an organic solvent to
mediate and engage in SEI chemistry. Rather, the salt alone must undergo the
reduction reactions that produce the SEI Steinrück  _et al._ (2020). In this
case, the specific chemical make-up of salt (specifically the anion), as well
as the precise positioning of species near the electrochemical interface will
strongly control the chemistry of SEI formation, as has been discussed at
length in literature Borodin  _et al._ (2017b); Steinrück  _et al._ (2018).
However, one aspect of SEI formation that is scarcely considered is the
thermodynamic activity of the salt ions that undergo this chemistry. The
activity, $a_{i}$, of the salt ions will strongly govern the driving force for
salt to undergo the SEI-forming chemical reactions. Our model allows us to
predict the thermodynamic activity of the salt as a function of the salt
concentration, the equations for which are given in the SI.

Figure 8: The difference in log of activity of water (blue) and salt (red)
from that at 0.2 m are plotted as a function of salt concentration for aqueous
LiTFSI electrolyte. Values of water activity are measured from molecular
dynamics explicitly (blue squares), and then fitted via a 4 th order
polynomial (blue dotted curve). The fitted curve is inserted into the Gibbs-
Duhem equation to reconstruct the activity of the salt (red dotted curve). The
kinks in the theoretical curves occur at the gel point. The theoretical curves
(smooth lines) were generated from Eqs. (S56), (S58) & (S59) in the SI, using
the parameters $\lambda=0.075$, $\xi_{+}=0.4$, $\xi_{-}=11.3$, $f_{+}=4$, and
$f_{-}=3$.

As we can see in Fig. 8, the activity of the LiTFSI salt
($a_{salt}=\sqrt{a_{+}a_{-}}$) is predicted by the theory (red smooth line) to
increase drastically with increasing salt concentration. In Fig. 8 we also
plot the activity of the LiTFSI salt as calculated from MD simulation. The
precise method of for computing the activity of the salt is outlined in the
SI. In brief, we explicitly measured the activity of water in the MD
simulation (blue squares) via the Osmotic Pressure Activity of Solvent (OPAS)
method Kohns  _et al._ (2016), which is explained in more detail in the SI.
Essentially, the activity of the water in an electrolyte of interest can be
determined from the osmotic pressure difference between a chamber with pure
water and the electrolyte of interest. We then fit a polynomial curve (blue
dotted line) to water activity as a function of salt concentration and
integrated the Gibbs-Duhem equation to reconstruct the activity of the salt
(red dotted line). We see that our reconstructed salt activity has a strong
qualitative agreement with the theoretical prediction. Note, the activities
plotted in Fig. 8 are plotted as a difference from a reference state
($\theta$) of 0.5 m. Of course, this data range does not deeply explore the
dilute regime Debye and Hückel (1923); Hückel (1925). For very dilute
electrolytes ($<$0.5 m), the ionic activity should decrease with
concentration. However, for the concentration range studied here, the salt
activity increases strongly as a function of salt concentration (at high salt
concentrations beyond the Debye-Hückel regime) primarily because Li+ starts to
become more associated to anions as opposed to solvent. Moreover, the
increased clustering reduces the overall translational degrees of freedom of
salt in the system. The theory also predicts a kink in chemical potential at
the gel point, leading to a more rapid increase in the salt chemical
potential, which is not observed in the reconstructed simulation curve. This
kink is likely due to the more rapid reduction of translational degrees of
freedom of the salt when increasing the salt concentration beyond the critical
gel concentration. Note, that we do not observe the this kink in our simulated
system. However, this could be due to the fitting process of the water
activity, which assumes a smooth transition through the gel point, and thus
yields a smooth curve for the Gibbs-Duhem reconstruction of the salt activity.
In order to resolve such a kink in the MD simulations, we would need a much
high density of simulation data near the gel point.

The overall increased salt activity will translate into an increase in the
reduction potential of any salt-derived SEI product. In order to demonstrate
this fact, consider the following reaction scheme, which was recently
determined to be the principle reaction pathway for LiF SEI formation in
LiTFSI WiSEs Steinrück  _et al._ (2020):

$\displaystyle\text{Li}^{+}+(\text{R-F})^{-}+\text{e}^{-}\rightleftharpoons\text{LiF}+R^{-}$
(16)

where for a LiTFSI electrolyte R=CF3-SO2-N-SO2-CF2. The equilibrium reduction
potential, $\Delta\Phi^{eq}$, is given by the Nernst equation

$\displaystyle\Delta\Phi^{eq}=E^{\theta}+k_{b}T\ln\left(\frac{a_{\text{Li}^{+}}a_{(\text{R-F})^{-}}}{a_{\text{LiF}}a_{\text{R}^{-}}}\right)$
(17)

where

$\displaystyle
E^{\theta}=\frac{\mu^{\theta}_{\text{Li}^{+}}\mu^{\theta}_{(\text{R-F})^{-}}}{\mu^{\theta}_{\text{LiF}}\mu^{\theta}_{\text{R}^{-}}}$
(18)

and $\mu^{\theta}_{i}$ is a reference chemical potential of species $i$
(assumed to be the standard dilute chemical potential for electrolytic
species). Because LiF is deposited on the electrode surface as a separate
phase, we can assume that $a_{\text{LiF}}$ is not _greatly_ affected by
electrolyte salt concentration and that $a_{\text{LiF}}\approx 1$. Our model
cannot definitively say much about the activity R-, but for simplicity we
approximate $a_{\text{R}^{-}}\approx 1$. In this case, the equilibrium
reduction potential for the LiF SEI formation will be shifted by the following
factor:

$\displaystyle\Delta\Phi^{eq}-E^{\theta}=k_{b}T\ln\left(a_{\text{Li}^{+}}a_{(\text{R-F})^{-}}\right)=2k_{b}T\ln
a_{salt}$ (19)

Thus, we can estimate that the increase in salt activity (from 0.2 m to 21 m)
would theoretically result in a positive shift of 0.69 eV (0.44 eV for MD)
toward higher electrode potentials, resulting in a strongly increased affinity
to form the SEI layer. In fact, this increase in affinity may explain why
purely-salt derived SEI layers require extreme salt concentrations: it is
simply thermodynamically unfavorable for the salt to be reduced at low salt
concentrations. Thus, understanding how the chemical potential of the salt
varies as a function of salt concentration may be a key consideration when
designing new WiSEs.

It has also been well-established that the coordination environment of Li+
strongly affects the SEI chemistry Xu  _et al._ (2007); Xu (2007); von Wald
Cresce  _et al._ (2012). For example, in Ref. 84, von Wald Cresce et al.
studied an organic ethylene carbonate (EC)/propylene carbonate (PC)
electrolyte linking the preference of PC over EC to be involved in SEI
formation to the preference of PC over EC to be in Li+ solvation sheath. In 21
m LiTFSI, the average coordination environment of Li+ contains roughly two
anions and two water molecules, as opposed to the typical pristine four-water
solvation shell in dilute solutions. Thus, the direct contact between Li+ and
the TFSI- is thought to enable their combined reduction to the LiF-rich SEI
layer Suo  _et al._ (2015). This is consistent with the thermodynamic
explanation based on salt activity presented above.

Furthermore, for a simple 3 component (cation, anion and solvent) electrolyte
the SEI chemistry that is possible is narrow, and can perhaps be grasped
intuitively without too much theoretical guidance. However, the next
generation of battery electrolytes may contain multiple solvents Wang  _et
al._ (2018); Zhang  _et al._ (2018); Dou  _et al._ (2018, 2019); Molinari and
Kozinsky (2020) or multiple salts Lui  _et al._ (2011); Suo  _et al._ (2016);
Kondou  _et al._ (2018); Molinari  _et al._ (2019a, b); Chen  _et al._
(2020b); Jiang  _et al._ (2020); Becker  _et al._ (2020). In this case,
predictive, chemical-specific models will be extremely important in building
our intuition. Knowing how the Li+ coordination environment changes for
different complex, super-concentrated salt/solvent blends will be extremely
non-trivial, and will help aid in understanding increasingly complex SEI
chemistry. Our current model, is equipped to handle the simple 3 component
WiSEs, but further development will be needed to generalize it to these more
complex mixtures.

Another aspect that can affect the electrochemical stability of WiSEs is the
activity of water, which factors strongly into the ability for the water to be
stabilized in the presence of large voltages. The activity of water is
expected to decrease with increasing salt concentration, as predicted in Fig.
8 by both our model (smooth line) and the simulation (smooth line) for LiTFSI
electrolyte. This is a result of water transitioning from being primarily
unbound to Li+ at low salt concentrations, to being primarily bound to Li+ at
high salt concentrations. The Li-water interaction is strongly energetically
favorable and reduces the activity of water. The Nernst equation again implies
that activity of water will tend to shift the oxidation potential of water by
a factor of $-4k_{B}T\log(a_{w})$, corresponding to a theoretical 0.37 eV
(0.21 eV from MD) shift in potential from a concentration of 0.2 m to a
concentration of 21 m. This implies that when the activity of water is reduced
below 1, the equilibrium potential for the oxygen evolution reaction is
shifted to more positive potentials. Additionally, it was also observed that
water molecules are considerably depleted from positive electrode surfaces
Vatamanu and Borodin (2017); McEldrew  _et al._ (2018) in super-concentrated
electrolytes. Thus, oxygen evolution reaction would be expected to be
suppressed even further than what was thermodynamically predicted here.

Of course electrochemical window may not be considered only in terms of
chemical thermodynamics. There are kinetic factors that will determine the
reactivity of water or the formation of SEI at active electrodes, such as e.g.
electron transfer reorganization energy, local populations of reactive
species, or electrocatalytic effects involving adsorbed intermediates.
Obviously such considerations lie beyond the scope of the presented analysis,
but even the ’classical’ concepts of activities discussed above, give a
qualitative indication of the trend in water reactivity and SEI formation with
the increase of the salt concentration, in line with available experimental
observations.

##  V Conclusion

In this article we have outlined the theory of ion aggregation and gelation in
water-in-salt electrolytes. Our model contains a small number of physically
transparent molecular parameters, which can be readily obtained from molecular
simulations. The model was validated across an extremely large concentration
range from dilute (0.2 m) to super-concentrated (21 m) for aqueous LiTFSI,
LiOTF, and LiFSI electrolytes. In our molecular simulations we observed
extensive ion clustering and the emergence of a percolating ion network, which
our model was able to represent as an ionic gel. The model is able to account
for the strong competition between ionic association and ion solvation, as was
observed in the performed simulations. In so doing, we are able to reproduce
the ion cluster distributions, Li+ coordination, and the degree of ion
aggregation as a function of salt concentration for all of the simulated
electrolytes.

Moreover, our model seems to be able to help explain various intriguing
properties of WiSE’s that have been observed in the literature. In particular,
there have been various observations detailing the heterogeneous nano-
structure in WiSEs, specifically the emergence of percolating ion networks
interpenetrated by water/Li+. Our theory is able to predict the conditions
required for such structures to emerge, as well as the general trend of
asymmetric clusters leading to cationic transference numbers increasing
considerably above 0.5 as salt concentration increases. Additionally, we were
able to use a thermodynamic argument, involving the simultaneous increase of
salt activity and decrease water activity at high salt concentrations, to
explain the preference of SEI formation over hydrogen evolution. This
preference of SEI formation is essential in extending the electrochemical
stability window of water-in-salt electrolytes. Thus, understanding this
process, may unlock the design principles for further optimization of WiSEs.

##  VI Acknowledgements

All authors would like to acknowledge the Imperial College-MIT seed fund. M.M.
and M.Z.B. acknowledge support from a Amar G. Bose Research Grant. Z.A.H.G was
supported through a studentship in the Centre for Doctoral Training on Theory
and Simulation of Materials at Imperial College London funded by the EPSRC
(EP/L015579/1) and from the Thomas Young Centre under grant number TYC-101.
A.A.K would like to acknowledge the research grant by the Leverhulme Trust
(RPG-2016- 223). This work used the Extreme Science and Engineering Discovery
Environment (XSEDE), which is supported by National Science Foundation grant
number ACI-1548562. S.B. was supported by the National Natural Science
Foundation of China (51876072) and the financial support from the China
Scholarship Council.

##  VII Appendix

###  VII.1 Mixture Stoichiometry

As described in the main chapter, we consider a polydisperse mixture of
$\sum_{lms}N_{lms}$ ionic clusters, each containing $l$ cations, $m$ anions,
and $s$ water molecules associated to cations (rank $lms$ cluster), and (if
present) an interpenetrating gel network containing $N_{+}^{gel}$ cations,
$N_{-}^{gel}$ anions, and $N_{0}^{gel}$ water molecules. We model the cations
to have a functionality (defined as the number of associations that the
species can make) of $f_{+}$, and anions to have a functionality of $f_{-}$.
This means that a cation is able to associate with $f_{+}$ anions or water
molecules. Likewise, anions may associate with up to $f_{-}$ cations. Water
molecules may associate (hydrate) to a single cation (functionality of one).

Our model treats the electrolyte as lattice fluid. We designate a single
lattice site to have the volume of a single water molecule, $v_{0}$. Thus, the
entire volume of the mixture, $V$, is divided into $\Omega=V/v_{0}$ lattice
sites. Moreover, cations will occupy $\xi_{+}=v_{+}/v_{0}$ lattice sites, and
anions will occupy $\xi_{-}=v_{-}/v_{0}$ lattice sites. Furthermore, when a
gel is formed, we distinguish between the volume fractions of gel (superscript
$gel$) and sol (superscript $sol$). The volume fractions in the sol and gel
constitutes the total volume fraction, $\phi_{j}$ of a given species, $j$, is
given by

$\displaystyle\phi_{j}=\phi_{j}^{sol}+\phi_{j}^{gel}$ (20)

in which the gel volume fractions is defined as
$\phi_{j}^{gel}=\xi_{j}N_{j}^{gel}/\Omega$, with $N_{j}^{gel}$ as the mole
number of species $j$ in the gel. The subscript $j=+,-,0$ corresponds to
cation, anion, and solvent, respectively. The sol volume fraction of cations,
anions, and solvent molecules have, respectively, the definitions

$\displaystyle\phi^{sol}_{+}=\sum_{lms}\xi_{+}lc_{lms}$ (21)
$\displaystyle\phi^{sol}_{-}=\sum_{lms}\xi_{-}mc_{lms}$ (22)
$\displaystyle\phi^{sol}_{0}=\sum_{lms}sc_{lms}$ (23)

where $c_{lms}=N_{lms}/\Omega$ is the dimensionless concentration of a $lms$
cluster (the number of $lms$ clusters per lattice site). Similarly, we define
$\phi_{\pm}=\phi_{+}+\phi_{-}$, which is the total volume fraction of the salt
in solution. For simplicity the mixture is assumed to be incompressible, i.e.

$\displaystyle 1=\phi_{\pm}+\phi_{0}=\phi_{+}+\phi_{-}+\phi_{0}$ (24)

$\phi_{+}$ and $\phi_{-}$ are not independent owing to electroneutrality:
$\phi_{+}/\xi_{+}=\phi_{-}/\xi_{-}$. The reduced volume of the mixture,
$\Omega$, can also be expressed in terms of the mole number of each
species/component due to the incompressibility constraint [Eq. (24)]

$\displaystyle\Omega=\sum_{lms}(\xi_{+}l+\xi_{-}m+s)N_{lmsq}+\xi_{+}N^{gel}_{+}+\xi_{-}N^{gel}_{-}+N^{gel}_{0}$
(25)

This definition must be used when differentiating the free energy of mixture.

###  VII.2 Pre-gel Cluster Distribution and Association Probabilities

Recall the the free energy of mixing, $\Delta F$, written in the main text:

$\displaystyle\beta\Delta F$
$\displaystyle=\sum_{l,m,s}\left[N_{lms}\log\left(\phi_{lms}\right)+N_{lms}\Delta_{lms}\right]$
$\displaystyle+\Delta^{gel}_{+}N^{gel}_{+}+\Delta^{gel}_{-}N^{gel}_{-}+\Delta^{gel}_{0}N^{gel}_{0}$
(26)

where $\beta=1/k_{B}T$ is the inverse thermal energy, and $\Delta_{lms}$ is
the free energy of formation for clusters of rank $lms$, and
$\Delta^{gel}_{i}$ is the free energy change of species $i$ upon association
to the gel.

We may differentiate the free energy with respect to $N_{lms}$ to obtain the
chemical potential of a cluster rank $lms$

$\displaystyle\beta\mu_{lms}$
$\displaystyle=\ln\phi_{lms}+1-(\xi_{+}l+\xi_{-}m+s)c_{tot}+\Delta_{lms}$
$\displaystyle+d\left[\phi_{0}(\xi_{+}l+\xi_{-}m+s)-s\right],$ (27)

where $c_{tot}=\sum_{lms}c_{lms}$ is the total reduced concentration, and
$d=\Delta^{gel\prime}_{+}c_{+}^{gel}+\Delta^{gel\prime}_{-}c_{-}^{gel}+\Delta^{gel\prime}_{0}c_{0}^{gel}$
(the $\prime$ notation refers to a derivative with respect to $\phi_{\pm}$),
and $c_{i}^{gel}=N_{i}^{gel}/\Omega$ is dimensionless concentration of $i$ in
the gel (number of species $i$ in the gel per lattice site). Additionally, we
may define the chemical potential of species immersed in the gel, which yields

$\displaystyle\beta\mu_{+}^{gel}=\Delta^{gel}_{+}-\xi_{+}c_{tot}+\xi_{+}d\phi_{0}$
(28)

for cations,

$\displaystyle\beta\mu_{-}^{gel}=\Delta^{gel}_{-}-\xi_{-}c_{tot}+\xi_{-}d\phi_{0}$
(29)

for anions, and

$\displaystyle\beta\mu_{0}^{gel}=\Delta^{gel}_{0}-c_{tot}+d\phi_{\pm}$ (30)

for water molecules.

The distribution of clusters can be derived by enforcing a chemical
equilibrium between all of the clusters and their bare constituents
(unassociated components)

$\displaystyle l[\text{bare cation}]+m[\text{bare anion}]+(s+q)[\text{bare
solvent}]\rightleftharpoons[lmsq\,\,\text{cluster}]$ (31)

Chemical equilibrium requires that the chemical potentials of bare species and
those in clusters are equivalent

$\displaystyle
l\mu_{100}+m\mu_{010}+s\mu_{001}=\mu_{lms}=l\mu^{+}_{lms}+m\mu^{-}_{lms}+s\mu^{0}_{lms}.$
(32)

In Eq. (32), we have defined the chemical potential of a cation, anion or
solvent molecule in an arbitrary cluster in the following manner

$\displaystyle\mu^{+}_{lms}=\frac{\partial\mu_{lms}}{\partial l}=\mu_{100}$
(33) $\displaystyle\mu^{-}_{lms}=\frac{\partial\mu_{lms}}{\partial
m}=\mu_{010}$ (34)

and

$\displaystyle\mu^{0}_{lms}=\frac{\partial\mu_{lms}}{\partial s}=\mu_{001}$
(35)

Thus, it becomes clear that the equilibrium condition in Eq. (32) requires
that species within an arbitrary cluster ($\mu^{i}_{lms}$, $i=+,-,0$) have
equivalent chemical potentials to the bare species not in the cluster
($\mu_{100}$, $\mu_{010}$, $\mu_{001}$).

We may solve Eq. (32) for an arbitrary cluster rank $lms$ obtaining the
following relationship:

$\displaystyle\phi_{lms}=K_{lms}\phi_{100}^{l}\phi_{010}^{m}\phi_{001}^{s}$
(36)

where $K_{lms}$ is the equilibrium constant with the following definition:

$\displaystyle K_{lms}=\exp(l+m+s-1-\Delta_{lms})$ (37)

As we mentioned in the main text, the number of water molecules in the cluster
($s$) is specified by the number of cations ($l$) and the number of anions
($m$) in the cluster: $s=f_{+}l-l-m+1$. Thus, strictly speaking, the
subscripts $lms$ describing the rank of a cluster for any variable (such as
$\Delta_{lms}$) could simply be described by the subscript $lm$. However, for
book-keeping purposed we will keep the subscript $lms$ throughout. The
partitioning of the species into clusters of different sizes is determined
almost entirely from $\Delta_{lms}$. As such, this is where much of the
physics of the ion/solvent association will be included. $\Delta_{lms}$
contains three major contributions: 1) _combinatorial_ describing the
multiplicity of clusters with the same number of constituents, 2) _binding_
describing the association enthalpy of the constituents in the cluster, and 3)
_configurational_ describing the configurational entropy change upon forming a
cluster from base constituents. As previously mentioned, in the sticky-cation
approximation we have neglected the non-ideal electrostatic contributions to
$\Delta_{lms}$, such as Debye-Hückle like screening or Born solvation free
energy. Thus,

$\displaystyle\Delta_{lms}=\Delta_{lms}^{comb}+\Delta_{lms}^{bind}+\Delta_{lms}^{conf}.$
(38)

The combinatorial contribution $\Delta^{comb}_{lms}$ is a purely entropic
contribution. Thus, it requires the combinatorial enumeration, $W_{lms}$, of
all of the possible ways a cluster with $l$ cations, $m$ anions, and $s$
solvent molecules can be formed:

$\displaystyle\Delta_{lms}^{comb}=-\log\left(W_{lms}\right)$ (39)

Here $W_{lms}$ can be derived in a two step procedure. First, we enumerate the
number of ways, $W_{lm}$ to construct a network containing $l$ anions and $m$
cations, which are associated together in an alternating fashion. This
combinatorial problem is well knownStockmayer (1952):

$\displaystyle
W_{lm}=\frac{(f_{+}l-l)!(f_{-}m-m)!}{l!m!(f_{+}l-l-m+1)!(f_{-}m-l-m+1)!}.$
(40)

In the second step, we enumerate the number of ways $s$ solvent molecules can
be placed on the cation-anion cluster. We know that we may only place the $s$
solvent molecules on the remaining $f_{+}l-l-m+1$ open cation sites, and
actually because there will be no open sites on the cations, and the solvent
molecules are indistinguishable, the cluster will contain $s=f_{+}l-l-m+1$
solvent molecules, with only one unique arrangement. Thus, $W_{lms}=W_{lm}$.

Next, the binding contribution, $\Delta^{bind}_{lms}$, can be describe simply
via the association free energies: $\Delta u_{+0}$ for cation-solvent
association and $\Delta u_{+-}$ for cation-anion association. For $l+m>0$, the
association free energy for an $lms$ cluster is

$\displaystyle\Delta^{bind}_{lms}=s\Delta u_{+0}+(l+m-1)\Delta u_{+-}$ (41)

The coefficient in front of $\Delta u_{+-}$ is due to the fact that there must
be that many cation/anion associations to form a cluster with $l$ cations and
$m$ anions. Again $s=f_{+}l-l-m+1$, thus

$\displaystyle\Delta^{bind}_{lms}=f_{+}l\Delta u_{+0}+(l+m-1)(\Delta
u_{+-}-\Delta u_{+0})$ (42)

For $l+m=0$, $s=1$ (free solvent molecule), we the binding energy is zero:

$\displaystyle\Delta^{bind}_{lms}=0.$ (43)

Recall, that our model does not allow for solvent molecules to form clusters
among themselves. For this reason, if a cluster contains 0 cations and anions,
the cluster will necessarily only contain a single solvent molecule,
corresponding to a free solvent molecule. Clearly, a free water molecule does
not form associations and thus $\Delta^{bind}_{001}=0$. Overall, we can write
$\Delta^{bind}_{lms}$ as

$\displaystyle\Delta^{bind}_{lms}=\left[f_{+}l\Delta u_{+0}+(l+m-1)(\Delta
u_{+-}-\Delta u_{+0})\right](1-\delta_{l,0}\delta_{m,0})$ (44)

where $\delta_{i,j}$ is the Kroenecker delta function.

Finally, for the configurational contribution, $\Delta^{conf}_{lms}$, as was
noted in ref. 42, we may model the associations as being semi-flexible. In
order to do so, associations are partitioned between trans and gauche
orientations, which differ by an energy (units $k_{B}T$) of $\epsilon$ (trans
conformation has the reference energy of 0, and gauche $\epsilon$). Flory’s
expression for the entropy of disorientation used in lattice fluid theoryFlory
(1942, 1953); Tanaka (1989); Tanaka and Ishida (1999) may be augmented to
account for the partitioning of associations into trans and gauche
orientations Flory (1956); Tanaka (2011):

$\displaystyle\Delta_{lm}^{conf}=-\ln\left(\frac{(\xi_{+}l+\xi_{-}m+s)\left[\frac{(Z-2)^{2g}}{Z2eg^{2g}(1-g)^{2(1-g)}}\right]^{l+m+s-1}}{\xi_{+}^{l}\xi_{-}^{m}}\right)+g\epsilon$
(45)

where $Z$ is the coordination number of the lattice, and $g$ is the fraction
of associations in the gauche conformation given byFlory (1956)

$\displaystyle g=\frac{(Z-2)e^{-\epsilon}}{1+(Z-2)e^{-\epsilon}}$ (46)

In the high temperature limit ($\epsilon\rightarrow 0$),
$g\rightarrow(Z-2)/(Z-1)$ for completely flexible chains. This yield’s Flory’s
original entropy of disorientation resultFlory (1942). In the low temperature
limit ($\epsilon\rightarrow\infty$), we obtain $g\rightarrow 0$, and all
associations are stuck in the trans orientation, thus associations become
completely inflexible. As we will see, when working within the Sticky-cation
approximation, the configurational entropy actually gets cancelled out, and
whether the associations are fully flexible, semi-flexible, or completely
inflexible does not matter.

Plugging in each contribution of $\Delta_{lms}$ into Eq. (36), we obtain the
distribution (for $l+m>1$)

$\displaystyle\tilde{c}_{lms}=\frac{W_{lm}}{\lambda_{+-}}\left(\frac{f_{+}\phi_{100}\lambda_{+-}}{\xi_{+}}\right)^{l}\left(\frac{f_{-}\phi_{010}\lambda_{+-}}{\xi_{-}}\right)^{m}\left(\phi_{001}\lambda_{+0}\right)^{f_{+}l-l-m+1}$
(47)

where $\lambda_{+-}$ and $\lambda_{+0}$ are the association constants for
cation-anion association and cation-solvent association, respectively. They
are defined generally as

$\displaystyle\lambda_{+-}=\frac{(Z-2)^{2g}}{Zg^{2g}(1-g)^{2(1-g)}}\exp\left(-\Delta
u_{ij}-g\epsilon\right)$ (48)

where $i=+$ and $j=-,0$ are indices denoting either cation-anion or cation-
solvent associations. For completely flexible associations,
$g\rightarrow(Z-2)/(Z-1)$ yielding

$\displaystyle\lambda^{flex}_{ij}=\frac{(Z-1)^{2}}{Z}\exp\left(-\Delta
u_{ij}\right).$ (49)

For completely inflexible associations, $g\rightarrow 0$ yielding

$\displaystyle\lambda^{inflex}_{ij}=\frac{1}{Z}\exp\left(-\Delta
u_{ij}\right).$ (50)

Note that in Eq. (47), we have removed the $s$, index, as the cluster is now
defined simply by the number of anions and cations (recall $s=f_{+}l-l-m+1$).
Furthermore, Eq. (47) gives the thermodynamically consistent number
distribution of clusters in terms of the volume fraction of the _bare_
cations, anions, and solvent molecules are known, which are generally unknown.

We may determine the bare species volume fractions by introducing the
association probabilities, $p_{ij}$, that an association site of species $i$
is occupied with an association to species $j$. These probabilities are useful
because we may write the reduced bare species’ volume fractions in terms of
these probabilities. For the bare cation volume fraction we have,

$\phi_{100}=\phi_{+}(1-p_{+-}-p_{+0})^{f_{+}}.$ (51)

The above equation arises because the probability that a given cation
association site will be ‘dangling’ (not participating in associations) will
be $1-p_{+-}-p_{+0}$. Thus for all $f_{+}$ sites to be dangling is
$(1-p_{+-}-p_{+0})^{f_{+}}$. Similarly, for the bare anion volume fraction we
have

$\phi_{010}=\phi_{-}(1-p_{-+})^{f_{-}},$ (52)

where $(1-p_{-+})^{f_{-}}$ is the probability that all $f_{-}$ anion sites
will be unassociated. Finally, for the bare solvent volume fraction we have

$\phi_{001}=\phi_{0}(1-p_{0+}),$ (53)

where $1-p_{0+}$ is the probability that a solvent molecule’s single
association site will not participate in an association. The probabilities are
related in the following manner due to conservation of associations:

$\displaystyle f_{+}\phi_{+}p_{+-}/\xi_{+}=f_{-}\phi_{-}p_{-+}/\xi_{-}=\zeta$
(54)

where $\zeta$ is the total number of anion-cation associations.

$\displaystyle f_{+}\phi_{+}p_{+0}/\xi_{+}=\phi_{0}p_{0+}=\Gamma$ (55)

where $\Gamma$ is the total number of cation-solvent associations. Tanaka also
employs the laws of mass action on the number of associations using the
association constants $\lambda_{+-}$ and $\lambda_{+0}$:

$\displaystyle\lambda_{+-}\zeta=\frac{p_{+-}p_{-+}}{(1-p_{+-}-p_{+0})(1-p_{-+})}.$
(56)

Similarly, for the cation-solvent associations we have

$\displaystyle\lambda_{+0}\Gamma=\frac{p_{+0}p_{0+}}{(1-p_{+-}-p_{+0})(1-p_{0+})}$
(57)

Here we are essentially treating $\lambda_{+-}$ and $\lambda_{+0}$ as
equilibrium constants for individual associations made. However, Eqs. (56) &
(57) are actually singular, because the “stickyness” of the cation dictates
that there will be no dangling sites on the cation:

$\displaystyle 1=p_{+-}+p_{+0}$ (58)

Therefore, we must divide Eqs. (56) & (57) to remove the singularity:

$\displaystyle\tilde{\lambda}=\frac{\lambda_{+-}}{\lambda_{+0}}=\frac{p_{-+}(1-p_{0+})}{p_{0+}(1-p_{-+})}.$
(59)

where $\tilde{\lambda}=\exp(-\Delta u_{+-}+\Delta)u_{+0}$, which is
independent of flexibility of associations discussed above.

The Eqs. (54)-(55) & Eqs. (58)-(59) provide a system from which we may solve
for $p_{+-}$, $p_{-+}$, $p_{0+}$, and $p_{+0}$ in terms of the overall species
volume fractions. The solutions are

$\displaystyle p_{+-}$
$\displaystyle=\frac{\psi_{0}-\phi_{+}+\tilde{\lambda}(\psi_{+}+\psi_{-})}{2(\tilde{\lambda}-1)\psi_{+}}$
$\displaystyle-\frac{\sqrt{4(\tilde{\lambda}-1)+(\tilde{\lambda}(\psi_{-}-\psi_{+})+\psi_{+}+\psi_{0})^{2}}}{2(\tilde{\lambda}-1)\psi_{+}},$
(60) $\displaystyle p_{-+}$
$\displaystyle=\frac{\psi_{0}-\phi_{+}+\tilde{\lambda}(\psi_{+}+\psi_{-})}{2(\tilde{\lambda}-1)\psi_{-}}$
$\displaystyle-\frac{\sqrt{4(\tilde{\lambda}-1)+(\tilde{\lambda}(\psi_{-}-\psi_{+})+\psi_{+}+\psi_{0})^{2}}}{2(\tilde{\lambda}-1)\psi_{-}},$
(61) $\displaystyle p_{+0}$
$\displaystyle=1-\frac{\psi_{0}-\phi_{+}+\tilde{\lambda}(\psi_{+}+\psi_{-})}{2(\tilde{\lambda}-1)\psi_{+}}$
$\displaystyle+\frac{\sqrt{4(\tilde{\lambda}-1)+(\tilde{\lambda}(\psi_{-}-\psi_{+})+\psi_{+}+\psi_{0})^{2}}}{2(\tilde{\lambda}-1)\psi_{+}},$
(62)

and

$\displaystyle p_{0+}$
$\displaystyle=\frac{\psi_{0}}{\psi_{+}}-\frac{\psi_{0}-\phi_{+}+\tilde{\lambda}(\psi_{+}+\psi_{-})}{2(\tilde{\lambda}-1)\psi_{0}}$
$\displaystyle+\frac{\sqrt{4(\tilde{\lambda}-1)+(\tilde{\lambda}(\psi_{-}-\psi_{+})+\psi_{+}+\psi_{0})^{2}}}{2(\tilde{\lambda}-1)\psi_{0}},$
(63)

where $\psi_{i}=f_{i}\phi_{i}/\xi_{i}=f_{i}\tilde{c}_{i}$. Note that
$\xi_{0}=f_{0}=1$.

Finally, we should rewrite the cluster distribution equation [Eq. (47)],
because the current cluster concentration distribution is written in terms of
the individual association constants ($\lambda_{+-}$ and $\lambda_{+0}$),
which will be singular and the “bare” cation volume fraction ($\phi_{100}$)
which will be zero. We may use Eqs. (51)-(59), to remove the singularities
from the $c_{lm}$ in Eq. (47). We obtain the following result:

$\displaystyle\tilde{c}_{lms}=\frac{\psi_{0}\alpha_{0}W_{lm}}{\tilde{\lambda}}\left(\tilde{\lambda}\frac{\psi_{+}\alpha_{+}}{\psi_{0}\alpha_{0}}\right)^{l}\left(\tilde{\lambda}\frac{\psi_{-}\alpha_{-}}{\psi_{0}\alpha_{0}}\right)^{m}(\tilde{\lambda}-\delta_{l,0}\delta_{m,0})$
(64)

where $\alpha_{0}=1-p_{0+}$, $\alpha_{+}=(1-p_{+-})^{f_{+}}$ and
$\alpha_{-}=(1-p_{-+})^{f_{-}}$ are the fraction of free water molecules,
cations and anions, respectively. Note that $\tilde{c}_{00}$ corresponds to
the dimensionless concentration of free water molecules. Equivalently, we may
write the cluster concentration distribution in terms of the ionic association
probabilities:

$\displaystyle\tilde{c}_{lm}=\tilde{c}_{salt}W_{lm}\mathcal{K}\left(\frac{p_{-+}}{(1-p_{-+})}(1-p_{+-})^{f_{+}-1}\right)^{l}\left(\frac{p_{+-}}{(1-p_{+-})}(1-p_{-+})^{f_{-}-1}\right)^{m}(\tilde{\lambda}-\delta_{l,0}\delta_{m,0})$
(65)

where $\tilde{c}_{salt}=\phi_{+}/\xi_{+}=\phi_{-}/\xi_{-}$ is the
dimensionless concentration of salt (# of salt molecules per lattice site),
and $\mathcal{K}=f_{\pm}(1-p_{-+})(1-p_{+-})/p_{\mp\pm}$. We may divide the
distribution by $2\tilde{c}_{salt}$ to obtain the cluster _probability_
distribution, $\alpha_{lm}$:

$\displaystyle\alpha_{lm}=W_{lm}\frac{\mathcal{K}}{2}\left(\frac{p_{-+}}{(1-p_{-+})}(1-p_{+-})^{f_{+}}\right)^{l}\left(\frac{p_{+-}}{(1-p_{+-})}(1-p_{-+})^{f_{-}}\right)^{m}(\tilde{\lambda}-\delta_{l,0}\delta_{m,0})$
(66)

which was similarly defined in Ref. 52.

We also note that the reduced association constant may now be written in the
following manner:

$\displaystyle\tilde{\lambda}=\exp\left(\Delta u_{+-}-\Delta
u_{+0}\right)=\exp\left(\beta(\Delta U_{+-}-\Delta U_{+0})\right)$ (67)

Thus, we can see that the entropic portion of the association constant
actually cancels out here, and is determined strictly by the energetics of
association. This, relies on the assumption that the entropy of solvent-cation
association is equivalent to cation-anion association. This assumption is
somewhat rough, and results from the Flory lattice configurational entropy. In
order to relax this assumption, we would have to modify Flory’s formulae. For
our purposes, however, this is unnecessary, as we see a satisfactory agreement
of the model with the simulations.

###  VII.3 More on the “Sticky” Cation Approximation

In order to probe the validity of the sticky cation approximation (SCA), we
may compare them to the general association probabilities derived from Eqs.
(54)-(57), which may be solved numerically. In Fig. 9, we show those
comparisons for different magnitudes of $\lambda_{+-}$ and $\lambda_{+0}$, but
with the same relative magnitude, thus keeping a constant $\tilde{\lambda}$ of
1/2. We see that when association constants are of order 102, the SCA matches
the numerical solutions almost exactly. Thus, we see clearly that this
approximation is valid for large association constants.

Figure 9: A comparison of the association probabilities using the full
numerical solution to eqs. (54)-(57) and the analytical association
probabilities (eqs. (60)-(63)) obtained via the “sticky” cation approximation
(SCA). In all panels, the association probabilities using the SCA correspond
to $\tilde{\lambda}=1/2$. The association constants, $\lambda_{+-}$ and
$\lambda_{+0}$, used to generate the numerical association probabilities are
given in the figures, and have an increasing magnitude going from the left to
right panels with the same relative ratio.

###  VII.4 Critical Gel Line

In the main text, we derived the expression for the critical gel line (the
critical gel concentration as a function of $\tilde{\lambda}$ or vice versa).
In principle, the association constant, $\tilde{\lambda}$, can be mapped to
the temperature of the electrolyte. However, it is important to note that when
we take a detailed look at the temperature dependence, we must also question
the applicability of the “Sticky-Cation Approximation” (SCA) across the
studied temperature range. As we mentioned in the previous section, in order
for the SCA to be appropriate, both $\lambda_{+-}$ and $\lambda_{+0}$ need to
be roughly 100 or more. However, we know that in general, association
constants will decrease exponentially as a function of temperature according
to the following relation:

$\displaystyle\lambda_{ij}=\exp\left(-\frac{\Delta
F_{ij}}{k_{b}T}\right)=\exp\left(-\frac{\Delta U_{ij}-T\Delta
S_{ij}}{k_{b}T}\right),$ (68)

where $\Delta F_{ij}$, $\Delta U_{ij}$, and $\Delta S_{ij}$ are respectively
the free energy, internal energy, and entropy of association between species
$i$ and $j$. Thus, it is clear that raising the temperature, could bring us
outside of conditions in which the SCA is appropriate. Furthermore, the
partitioning of the free energy of association between energetic and entropic
contributions will also be critical in determining the temperature dependence
of the association constants. Upon inspection of eqs. (48), (49), and (50) in
comparison to eq. (68), we see that the entropy of association can take on the
following forms depending on the flexibility of the associations. For semi-
flexible associations, we have

$\displaystyle\Delta
S_{ij}=k_{B}\ln\left(\frac{(Z-2)^{2g}}{Zg^{2g}(1-g)^{2(1-g)}}\right).$ (69)

For completely flexible associations, we have

$\displaystyle\Delta S^{flex}_{ij}=k_{B}\ln\left(\frac{(Z-1)^{2}}{Z}\right).$
(70)

Finally, for completely inflexible associations, we have

$\displaystyle\Delta S^{inflex}_{ij}=-k_{B}\ln\left(Z\right)$ (71)

Thus, for fully flexible, associations the entropy of association is strictly
positive, and for fully inflexible associations the entropy of association is
strictly negative. As we have explained, within the SCA, the entropy of
association will always cancel out, and it will not affect the temperature
dependence of $\tilde{\lambda}$ or the critical gel line. However, if the SCA
is relaxed, the magnitude of the entropy of association will heavily affect
the temperature dependence of $\lambda_{ij}$, and thus will affect the
critical gel line.

In order to elucidate this point, we can compute the gel line in the
temperature-concentration plane for LiTFSI electrolyte. In the main text,
$\tilde{\lambda}$ was found to be 0.075 for LiTFSI electrolyte at 300K. We may
also compute the values of the individual association constants using
association probabilities computed from MD simulation (from eqs. (56) and
(57)). For LiTFSI electrolyte at 300K, we found that $\lambda_{+-}=159$ and
$\lambda_{+0}=2052$. Thus, as was explained in the previous section, these
large association constants clearly indicate that the SCA will be appropriate
for our system. However, when the temperature is increased both $\lambda_{+-}$
and $\lambda_{+0}$ decrease depending on the magnitude of the entropy of
association.

Figure 10: a) We plot the association constants for cation–anion association
($\lambda_{+-}$, solid, colored curves) and cation–solvent association
($\lambda_{+0}$, dotted, colored curves) as a function of temperature, for a
range of different association entropies (2-4$k_{B}$). Additionally, we plot
the reduced association constant ($\tilde{\lambda}$, black curve, used in the
sticky cation approximation) as a function of temperature. b) We plot the
critical gel line in the temperature–concentration plane, using the sticky
cation approximation (black curve), and the general case (allowing for empty
association sites on the cation) for a range of different association
entropies (2-4$k_{B}$).

This can be seen in fig. 10a, where the $\lambda_{+-}$ (colored solid curves)
and $\lambda_{+0}$ (colored dotted curves) are plotted as a function of
temperature for various association entropies ranging from $-4k_{B}$–$2k_{B}$.
Generally, it can be seen in fig. 10a that as the entropy of association
decreases, the the association constants will decreases more drastically as a
function of increasing temperature. For reference, in ref. McEldrew  _et al._
(2018), the temperature dependent studies of clustering in emimBF4 yielded an
association entropy of -3.3$k_{B}$. A full temperature dependent study of
clustering would be necessary to determine the entropy of association for the
WiSEs studied here. On the contrary, $\tilde{\lambda}$ actually increases as a
function of temperature–in the high temperature limit $\tilde{\lambda}$ would
tend toward 1.

For increasing temperature, the SCA actually predicts increasing ionic
association. Thus, the black curve in fig. 10b shows that the critical gel
concentration mildly decreases as a function of increasing temperature. This
is a counter-intuitive result. Essentially, because the SCA conserves
associations (e.g. every cation–anion association that is formed requires a
cation–solvent association to be broken) the overall translational degrees of
freedom of the mixture are conserved whether the electrolyte favors
cation–anion or cation–solvent association. This means that the SCA predicts
that there will be no entropic driving force to favor either type of
association. Translational entropy is not increased by forming or breaking ion
associations. Instead, for the SCA, increasing temperature only scales the
energetic driving force for associations. Because $\Delta U_{+0}<\Delta
U_{+-}$¡0, the reduced cation–solvent association energy ($\Delta
u_{+0}=\beta\Delta U_{+0}$) increases more rapidly as a function of
temperature than the reduced cation–anion association energy ($\Delta
u_{+-}=\beta\Delta U_{+-}$).

On the other hand, if we relax the SCA, we see that association (both cation-
water and cation-anion association) decreases as a function of increasing
temperature. This results in the critical gel concentration increasing with
temperature, when the entropy of association is negative (expected for
inflexible associations).

###  VII.5 Post Gel Association Probabilities

Once the critical gel threshold has been surpassed, we must determine the
volume fractions of cations, anions, and water in the gel, $\phi_{+}^{gel}$,
$\phi_{-}^{gel}$, and $\phi_{0}^{gel}$ respectively. In order to do this, we
employ Flory’s treatment of the post-gel regime in which the volume fraction
of free ions can be written equivalently in terms of overall association
probabilities, $p_{ij}$, and association probabilities taking into account
only the species residing in the sol, $p^{sol}_{ij}$

$\displaystyle\phi_{+}(1-p_{+-})^{f_{+}}=\phi_{+}^{sol}(1-p^{sol}_{+-})^{f_{+}}$
(72)
$\displaystyle\phi_{-}(1-p_{-+})^{f_{-}}=\phi_{-}^{sol}(1-p^{sol}_{-+})^{f_{-}}$
(73) $\displaystyle\phi_{0}(1-p_{0+})=\phi_{0}^{sol}(1-p^{sol}_{0+})$ (74)

Where $\phi_{i}^{sol}=1-\phi_{i}^{gel}$ is the volume fraction of species $i$
in the sol. We may determine each of the three unknown sol volume fractions
($\phi_{+}^{sol}$, $\phi_{-}^{sol}$, $\phi_{0}^{sol}$), as well as the four
unknown sol association probabilities ($p_{+-}^{sol}$, $p_{-+}^{sol}$,
$p_{+0}^{sol}$, $p_{0+}^{sol}$) using Eqs. (72)-(74) in addition to Eqs.
(54)-(55) & Eqs. (58)-(59), using sol-specific quantities in this case.

Thus, we have 7 equations and 7 unknowns (four sol association probabilities
and three sol species volume fractions). The fraction of species, $i$, in the
gel is simply given by

$\displaystyle w_{i}^{gel}=1-\phi^{sol}_{i}/\phi_{i}$ (75)

Note that prior to the critical gel concentration, we have the trivial
solution that $p_{\pm\mp}=p^{sol}_{\pm\mp}$ and $\phi_{i}=\phi^{sol}_{i}$,
yielding a gel fraction of $w_{i}^{gel}=0$. However, beyond the gel point,
there is a non-trivial solution yielding $w_{i}^{gel}>0$.

###  VII.6 Theoretical Expressions for Activity

Before deriving our model’s expressions for activity it is useful to define
the various related thermodynamic quantities. In general, the chemical
potential of any species $i$ can be parsed as the following:

$\displaystyle\mu_{i}$ $\displaystyle=\mu_{i}^{\theta}+k_{B}T\ln a_{i}$
$\displaystyle=\mu_{i}^{\theta}+k_{B}T\ln\left(\frac{c_{i}}{c_{i}^{\theta}}\right)+\mu_{i}^{ex}$
$\displaystyle=\mu_{i}^{\theta}+k_{B}T\ln\left(\frac{c_{i}}{c_{i}^{\theta}}\right)+k_{B}T\ln\gamma_{i}$
(76)

where $\mu_{i}^{\theta}$ is the reference state chemical potential, $a_{i}$ is
the activity, $\mu_{i}^{ex}$ is the excess chemical potential, $\gamma_{i}$ is
the activity coefficient, $c_{i}$ is the concentration (arbitrary units) and
$c_{i}^{\theta}$ is the concentration in the reference state of species $i$.
The activity is thus broken down as $a_{i}=\gamma_{i}c_{i}/c_{i}^{\theta}$,
and the excess chemical potential has the definition
$\mu_{i}^{ex}=k_{B}T\ln\gamma_{i}$. For ionic species this reference state
chemical potential is taken to be at “infinite dilution”. Of course at
infinite dilution we run into the problem that the concentration is zero, and
thus the chemical potential would tend toward $-\infty$. In order to avoid
this, it is standard to define the reference concentration to be 1 (arbitrary
units). In our case, the reference chemical potentials are defined prior to
mixing, and thus the reference concentration corresponds to species’ volume
fractions of 1. The chemical potential of a rank $lms$ cluster was written in
Eq. (27). As in Ref. 52, we may write the ionic activities using Eq. (27) and
subtract off non-zero excess chemical potential

$\displaystyle\log
a_{+}=\beta(\mu_{100}-\mu^{\theta}_{100})=\beta\mu_{100}-\ln\left\\{(1-p^{\theta}_{+-}-p^{\theta}_{+0})^{f_{+}}\right\\}-(1-\xi_{+}c^{\theta}_{tot})-\xi_{+}\phi^{\theta}_{0}d^{\theta}$
(77)

for cations and

$\displaystyle\log
a_{-}=\beta(\mu_{010}-\mu^{\theta}_{010})=\beta\mu_{010}-\ln\left\\{(1-p^{\theta}_{-+})^{f_{-}}\right\\}-(1-\xi_{-}c^{\theta}_{tot})-\xi_{-}\phi^{\theta}_{0}d^{\theta}$
(78)

for anions where the $\theta$ again corresponds to the value at infinite
dilution, $p^{\theta}_{-+}=p^{\theta}_{+-}=d^{\theta}=0$ and
$\phi_{0}^{\theta}=c^{\theta}_{tot}=1$. Thus, we may write the ionic
activities explictly as:

$\displaystyle\ln
a_{+}=\ln\phi_{+}+f_{+}\ln\left\\{\frac{(1-p_{+-}-p_{+0})}{1-p^{\theta}_{+-}-p^{\theta}_{+0}}\right\\}+\xi_{+}(1-c_{tot}+\phi_{0}d)$
(79)

for cations and

$\displaystyle\ln
a_{-}=\ln\phi_{-}+f_{-}\ln\left\\{(1-p_{-+})\right\\}+\xi_{-}(1-c_{tot}+\phi_{0}d)$
(80)

for anions. In Eq. (79), we have the ratio of two divergent terms:
$(1-p_{+-}-p_{+0})/(1-p^{\theta}_{+-}-p^{\theta}_{+0})$. This ratio is finite,
however, and can be specified via Eq. (57) as the following:

$\displaystyle\frac{(1-p_{+-}-p_{+0})}{1-p^{\theta}_{+-}-p^{\theta}_{+0}}=\frac{p_{+0}}{1-p_{0+}}=\frac{1-p_{+-}}{1-p_{0+}}$
(81)

Thus, the cation activity is written as

$\displaystyle\ln
a_{+}=\ln\phi_{+}+f_{+}\ln\left\\{\frac{1-p_{+-}}{1-p_{0+}}\right\\}+\xi_{+}(1-c_{tot}+\phi_{0}d)$
(82)

The activity of water, $a_{0}$, is simply expressed as

$\displaystyle\ln
a_{0}=\beta\mu_{001}=\ln\phi_{0}+\ln\left\\{(1-p_{+0})\right\\}+1-c_{tot}+\phi_{\pm}d.$
(83)

In Eq. (83), we do not need to adjust the activity to the infinite dilution
reference state as with the ions, because the mixing reference state for water
is already at the appropriate infinite dilution reference state. The final
thing we must specify is $d$, which we recall has the definition:
$d=\Delta^{gel\prime}_{+}c_{+}^{gel}+\Delta^{gel\prime}_{-}c_{-}^{gel}+\Delta^{gel\prime}_{0}c_{0}^{gel}$,
where the $\prime$ notiation corresponds to a derivative with respect to
$\phi_{\pm}$. Thus we need to know the functional dependence of
$\Delta^{gel}_{i}$. This is accomplished by enforcing an equilibrium between
species in the sol and species in the gel ($\mu_{i}=\mu_{i}^{gel}$) and
solving for $\Delta_{i}^{gel}$ (again making the chemical potential reference
state that at infinite dilution):

$\displaystyle\Delta^{gel}_{+}(\phi_{\pm})=\ln\phi_{+}+f_{+}\ln\left\\{\frac{1-p_{+-}}{1-p_{0+}}\right\\}+\xi_{+}$
(84)
$\displaystyle\Delta^{gel}_{-}(\phi_{\pm})=\ln\phi_{-}+f_{-}\ln\left\\{(1-p_{-+})\right\\}+\xi_{-}$
(85) $\displaystyle\Delta^{gel}_{0}(\phi_{\pm})$
$\displaystyle=\ln\phi_{0}+\ln\left\\{(1-p_{0+})\right\\}+1$ (86)

Thus, we may simply plug in these expressions to $d$, obtaining the
thermodynamic activities of each species used to generate the theoretical
curves in Fig. 8 in the main text.

###  VII.7 Additional Results for LiOTF and LiFSI

In the main text, we discussed molecular simulation results and comparisons
with the developed theory, but we primarily focused on LiTFSI electrolyte.
Here we will outline additional results and comparisons for LiOTF and LiFSI
electrolyte, as well.

Figure 11: The cluster size probability distributions are plotted for a
manifold salt concentrations for (a) LiOTF and (c) LiFSI electrolytes. The
cluster size cumulative distributions are also plotted for a manifold salt
concentrations for (b) LiOTF and (d) LiFSI electrolytes.

In Fig. 11, we show the probability and cumulative distribution of clusters as
a function of cluster size (total number of ions) for various salt
concentrations for both LiOTF (Figs. 11a&b) and LiFSI (Figs. 11c&d). The
general trends for the LiOTF and LiOTF electrolytes are very similar to the
LiTFSI electrolyte discussed in the main text. In both LiFSI and LiTFSI
electrolyte, there is a clear secondary peak in the cluster probability
distribution indicating the presence of a percolating ion network when the
salt concentration is 21 m. For LiFSI, there also is a subtle secondary peak
(located at a cluster size of roughly $10^{2}$ ions) for 15 m salt
concentration. As we have mentioned in the main text, such a peak indicates
the existance of a percolating cluster for 15 m LiFSI, which is not expected
from the theory until about 21m, as LiFSI was found to be the least
associating salt that we have studied. A possible explanation for the
observation of percolation at such a low concentration in LiFSI could be the
finite simulation box size, which would tend to induce percolation prematurely
as compared to an infinite or macroscopic system.

Figure 12: The probability distributions of cluster sizes (total number of
ions in the cluster) are plotted for various salt concentrations in (a) LiOTF
and (c) LiFSI electrolyte. The cumulative distributions of cluster sizes
(total number of ions in the cluster) are plotted for various salt
concentrations in (b) LiOTF and (d) LiFSI electrolyte. Theoretical parameters
for LiOTF were $\lambda=0.184$, $\xi_{+}=0.4$, $\xi_{-}=5.9$, $f_{+}=4$, and
$f_{-}=3$. Theoretical parameters for LiFSI were $\lambda=0.049$,
$\xi_{+}=0.4$, $\xi_{-}=6.8$, $f_{+}=4$, and $f_{-}=3$

In Fig. 12, we show theory and MD simulation comparisons for LiOTF electrolyte
(Fig. 12a-c) and for LiFSI electrolyte (Fig. 12d-f). Overall, the theory is
able to quantitatively reproduce the ion cluster fractions (free ions, ion
pairs, high order finite clusters, and gel), the Li+ coordination numbers, and
degree of aggregation as functions of salt concentration for both LiOTF and
LiFSI electrolytes. Recall, our theoretical curves required the computation of
the association constant $\tilde{\lambda}$ from simulations via the ionic
association probabilities plotted in Fig. 4 in the main text.

###  VII.8 Ion Association and Free Ion Lifetimes

In Fig. 13, we plot the mean lifetimes of ion associations and free ions as a
function of salt concentration for aqueous LiTFSI, LiFSI, and LiOTF
electrolytes.

Figure 13: The mean lifetimes in picoseconds of (a) ion associations and (b)
free ions are plotted as a function of salt concentration for aqueous LiTFSI
(blue), LiFSI (red), and LiOTF (yellow) electrolytes

###  VII.9 Molecular Dynamics Simulation Methodology

In this study, we performed all-atom classical MD simulations using LAMMPS
Plimpton (1995). We performed a set of simulations for three different salts:
LiTFSI, LiFSI and LiOTF. For LiFSI and LiOTF salts we performed a series of
fully periodic simulations at molal concentrations of 0.28 m, 0.5 m, 1 m, 2 m,
4 m, 7 m, 10 m, 12 m, 15 m, and 21 m. For LiTFSI salt we peformed the periodic
simulations at molal concentrations of 0.28 m, 0.5 m, 1 m, 2 m, 4 m, 7 m, 10
m, 12 m, 15 m, 16 m, 17 m, 18 m, 19 m, 20 m, and 21 m.

Each simulation contained 1000 water molecules and enough ion pairs to created
as close to the specifies concentrations (0.2 m-21 m) as possible. The
simulations were performed at fixed temperature (300 K) and pressure (1 bar),
with Nose-Hoover thermostat and barostat until the density of the fluid
equilibrated (10 ns with 1 fs time steps) and then, production runs were
performed for 20 ns. The initial configurations for all simulations were
generated using the open-source software, PACKMOL Martínez  _et al._ (2009).
All MD simulations were visualized using the open-source software, VMD
Humphrey  _et al._ (1996).

For all ionic species we employed the CL$\&$P force field, which was developed
for ionic liquid simulations, with same functional form as the OPLSAA force
field Lopes and Pádua (2012). Given the dense ionic nature of our systems, we
expect the CL$\&$P force field to be appropriate for WiSEs. For water, we
employed the spc/e force field. Interatomic interactions are determined using
Lorentz-Berthelot mixing rules. Long range electrostatic interactions were
computed using the Particle-Particle Particle-Mesh (PPPM) solver (with a cut-
off length of 12 $\AA$).

###  VII.10 Computing Activity from Molecular Dynamics

In order to measure the activity of the species in our MD simulations, we
employed the OPAS (osmotic pressure for activity of solvents) method developed
in Ref. 81, which involves measuring the osmotic pressure difference between a
chamber containing only solvent molecules, and a chamber of an electrolyte of
interest (see Fig. 14). The electrolyte is described by the same force field
as in the previously discussed fully periodic simulations. However, the
chambers are separated by an FCC sheet of hard spheres, which prevent salt
molecules from permeating, but allow solvent molecules to pass freely in
between the chambers. The chamber dimensions are given in the caption of Fig.
14. The system pictured in Fig. 14, is fully periodic, but the dynamics is
performed in an NVT ensemble, with a Nose-Hoover thermostat at a target
temperature of 300 K. The system is simulated until the internal pressure and
concentration distributions have equilibrated, and we designated 10 ns for
this equilibration for each simulation. We then perform production runs of 10
ns for each simulated concentration in which the average pressure in each
chamber is calculated. The initial number of water molecules in the water-only
chamber is 1000 for all simulations. The number of ions and water molecules in
the electrolyte chamber are given in Tbl. 1 for each simulated concentration.

Conc. (m) |  |  $\\#$ of ion pairs |  |  $\\#$ of waters
---|---|---|---|---
0.5m |  | 9 |  | 1187
1m |  | 18 |  | 1117
2m |  | 36 |  | 1000
4m |  | 72 |  | 822
7m |  | 81 |  | 642
10m |  | 97 |  | 539
12m |  | 105 |  | 485
15m |  | 113 |  | 418
21m |  | 125 |  | 330
Table 1: The number of ion pairs and water molecules initially placed in the
electrolyte chamber for the OPAS method simulations.

The difference in average pressure yields the osmotic pressure $\Pi$, which
can be related directly to the activity of the water, $a_{0}$ in the
electrolyte chamber (using the Gibbs-Duhem Relation for the pure water chamber
and an assumption of chemical equilibrium between the two chambers):

$\displaystyle\log(a_{0})=-\frac{\Pi}{\rho_{0}k_{B}T}$ (87)

where $\rho_{0}$ is the molar concentration of water in the chamber containing
only water.

Figure 14: A schematic of the MD simulation box set-up used to employ the OPAS
method for computing solvent activity. The dimensions of the chamber on the
right containing solely water are 30 $\AA$ $\times$ 30 $\AA$ $\times$ 33.2
$\AA$. The box on the left containing the electrolyte has the dimensions 30
$\AA$ $\times$ 30 $\AA$ $\times$ 41.7 $\AA$.

Upon computing the activity of the water in the electrolyte as a function of
salt concentration, we can use the Gibbs-Duhem relation for the electrolyte
chamber to compute chemical potential of the salt (with respect to a dilute
reference state):

$\displaystyle\Delta\mu_{s}(x)=-\frac{1}{2}\int\frac{1-x}{x}d\mu_{0}(x)$ (88)

where $x$ is the mole fraction of salt, and $\mu_{s}=(\mu_{+}+\mu_{-})/2$ is
the mean chemical potential of the salt. $\Delta\mu_{s}(x)$ is the difference
mean chemical potential of the salt from a reference salt concentration of
0.5m and that at a mole fraction of $x$. We can numerically compute the
integral in Eq. (88) using the simulated values of the activity of water
[$d\mu_{0}=d\ln(a_{w})$], however our numerical accuracy is greatly limited by
the discreteness of the activity measurements. For this reason, we instead fit
a polynomial (order 4) to the simulated activity and compute the integral
exactly, yielding the “MD Gibbs-Duhem” curve in Fig. 8 in the Main text.

## References

  * Suo  _et al._ (2013) L. Suo, Y.-S. Hu, H. Li, M. Armand, and L. Chen, Nat. Commun. 4, 1481 (2013).
  * Sodeyama  _et al._ (2014) K. Sodeyama, Y. Yamada, K. Aikawa, A. Yamada, and Y. Tateyama, J. Phys. Chem. C 118, 14091 (2014).
  * Suo  _et al._ (2015) L. Suo, O. Borodin, T. Gao, M. Olguin, J. Ho, X. Fan, C. Luo, C. Wang, and K. Xu, Science 350, 938 (2015).
  * Yamada  _et al._ (2016) Y. Yamada, K. Usui, K. Sodeyama, S. Ko, Y. Tateyama, and A. Yamada, Nat. Energy 1, 16129 (2016).
  * Wang  _et al._ (2016) J. Wang, Y. Yamada, K. Sodeyama, C. H. Chiang, Y. Tateyama, and A. Yamada, Nat. Commun. 7, 12032 (2016).
  * Gambou-Bosca and Bélanger (2016) A. Gambou-Bosca and D. Bélanger, Journal of Power Sources 326 (2016), 10.1016/j.jpowsour.2016.04.088.
  * Sun  _et al._ (2017) W. Sun, L. Suo, F. Wang, N. Eidson, C. Yang, F. Han, Z. Ma, T. Gao, M. Zhu, and C. Wang, Electrochem. Commun. 82, 71 (2017).
  * Suo  _et al._ (2017) L. Suo, O. Borodin, Y. Wang, X. Rong, W. Sun, X. Fan, S. Xu, M. A. Schroeder, A. V. Cresce, F. Wang, C. Wang, Y.-S. Hu, K. Xu, and C. Wang, Advanced Energy Materials 7, 1701189 (2017).
  * Dong  _et al._ (2017a) X. Dong, H. Yu, Y. Ma, J. L. Bao, D. G. Truhlar, Y. Wang, and Y. Xia, Chemistry - A European Journal 23 (2017a), 10.1002/chem.201700063.
  * Diederichsen  _et al._ (2017) K. M. Diederichsen, E. J. McShane, and B. D. McCloskey, ACS Energy Letters 2 (2017), 10.1021/acsenergylett.7b00792.
  * Yang  _et al._ (2017a) C. Yang, J. Chen, T. Qing, X. Fan, W. Sun, A. von Cresce, M. S. Ding, O. Borodin, J. Vatamanu, M. A. Schroeder, N. Edison, C. Wang, and K. Xu, Joule 1, 122 (2017a).
  * Yang  _et al._ (2017b) C. Yang, L. Suo, O. Borodin, F. Wang, W. Sun, T. Gao, X. Fan, S. Hou, Z. Ma, K. Amine, K. Xu, and C. Wang, Proceedings of the National Academy of Sciences of the United States of America 114 (2017b), 10.1073/pnas.1703937114.
  * Wang  _et al._ (2018) F. Wang, O. Borodin, M. S. Ding, M. Gobet, J. Vatamanu, X. Fan, T. Gao, N. Edison, Y. Liang, W. Sun, S. Greenbaum, K. Xu, and C. Wang, Joule 2, 927 (2018).
  * Leonard  _et al._ (2018) D. P. Leonard, Z. Wei, G. Chen, F. Du, and X. Ji, ACS Energy Lett. 3, 373 (2018).
  * Yamada  _et al._ (2019) Y. Yamada, J. Wang, S. Ko, E. Watanabe, and A. Yamada, Nature Energy 4, 269 (2019).
  * Yang  _et al._ (2019) C. Yang, J. Chen, X. Ji, T. P. Pollard, X. Lu, C. J. Sun, S. Hou, Q. Liu, C. Liu, T. Qing, Y. Wang, O. Borodin, Y. Ren, K. Xu, and C. Wang, “Aqueous Li-ion battery enabled by halogen conversion–intercalation chemistry in graphite,”  (2019).
  * Dou  _et al._ (2019) Q. Dou, Y. Lu, L. Su, X. Zhang, S. Lei, X. Bu, L. Liu, D. Xiao, J. Chen, S. Shi, and X. Yan, Energy Storage Materials 23 (2019), 10.1016/j.ensm.2019.03.016.
  * Lewis  _et al._ (2020) N. H. Lewis, Y. Zhang, B. Dereka, E. V. Carino, E. J. Maginn, and A. Tokmakoff, The Journal of Physical Chemistry C (2020).
  * Borodin  _et al._ (2017a) O. Borodin, L. Suo, M. Gobet, X. Ren, F. Wang, A. Faraone, J. Peng, M. Olguin, M. Schroeder, M. S. Ding, E. Gobrogge, A. von Wald Cresce, S. Munoz, J. A. Dura, S. Greenbaum, C. Wang, and K. Xu, ACS nano 11, 10462 (2017a).
  * Borodin  _et al._ (2020) O. Borodin, J. Self, K. A. Persson, C. Wang, and K. Xu, Joule 4, 69 (2020).
  * Li  _et al._ (2020) M. Li, C. Wang, Z. Chen, K. Xu, and J. Lu, Chemical Reviews (2020).
  * Chen  _et al._ (2020a) M. Chen, G. Feng, and R. Qiao, Current Opinion in Colloid & Interface Science 47, 99 (2020a).
  * Vatamanu and Borodin (2017) J. Vatamanu and O. Borodin, J. Phys. Chem. Lett. 8, 4362 (2017).
  * McEldrew  _et al._ (2018) M. McEldrew, Z. A. Goodwin, A. A. Kornyshev, and M. Z. Bazant, The journal of physical chemistry letters 9, 5840 (2018).
  * Li  _et al._ (2018) Z. Li, G. Jeanmairet, T. Méndez-Morales, B. Rotenberg, and M. Salanne, The Journal of Physical Chemistry C 122, 23917 (2018).
  * Choi  _et al._ (2018a) J.-H. Choi, H. Lee, H. R. Choi, and M. Cho, Annual review of physical chemistry 69, 125 (2018a).
  * Lim  _et al._ (2018) J. Lim, K. Park, H. Lee, J. Kim, K. Kwak, and M. Cho, Journal of the American Chemical Society 140, 15661 (2018).
  * Yu  _et al._ (2020) Z. Yu, L. A. Curtiss, R. E. Winans, Y. Zhang, T. Li, and L. Cheng, The Journal of Physical Chemistry Letters 11, 1276 (2020).
  * Andersson  _et al._ (2020) R. Andersson, F. Årén, A. A. Franco, and P. Johansson, Journal of the Electrochemical Society 167, 140537 (2020).
  * Debye and Hückel (1923) P. Debye and E. Hückel, Physikalische Zeitschrift 24, 185 (1923).
  * Hückel (1925) E. Hückel, Phys. Z 26, 93 (1925).
  * Fedorov and Kornyshev (2014) M. V. Fedorov and A. A. Kornyshev, Chem. Rev. 114, 2978 (2014).
  * Dong  _et al._ (2017b) K. Dong, X. Liu, H. Dong, X. Zhang, and S. Zhang, Chem. Rev. 117, 6636 (2017b).
  * Fawcett and Tikanen (1996) W. R. Fawcett and A. C. Tikanen, The Journal of Physical Chemistry 100, 4251 (1996).
  * Abbas  _et al._ (2007) Z. Abbas, E. Ahlberg, and S. Nordholm, Fluid phase equilibria 260, 233 (2007).
  * Vincze  _et al._ (2010) J. Vincze, M. Valiskó, and D. Boda, The Journal of chemical physics 133, 154507 (2010).
  * Shilov and Lyashchenko (2015) I. Y. Shilov and A. K. Lyashchenko, The Journal of Physical Chemistry B 119, 10087 (2015).
  * Bjerrum (1926) N. Bjerrum, “K. danske vidensk. selsk.” (1926).
  * Ebeling and Grigo (1980) W. Ebeling and M. Grigo, Annalen der Physik 492, 21 (1980).
  * Levin and Fisher (1996) Y. Levin and M. E. Fisher, Physica A: Statistical Mechanics and its Applications 225, 164 (1996).
  * Marcus and Hefter (2006) Y. Marcus and G. Hefter, Chemical reviews 106, 4585 (2006).
  * McEldrew  _et al._ (2020a) M. McEldrew, Z. A. H. Goodwin, H. Zhao, M. Z. Bazant, and A. A. Kornyshev, arXiv:2010.02091 (2020a).
  * Chen  _et al._ (2018) M. Chen, Z. A. H. Goodwin, G. Feng, and A. A. Kornyshev, J. Electroanal. Chem. 819, 347 (2018).
  * Feng  _et al._ (2019) G. Feng, M. Chen, S. Bi, Z. A. Goodwin, E. B. Postnikov, N. Brilliantov, M. Urbakh, and A. A. Kornyshev, Physical Review X 9, 021024 (2019).
  * Goodwin and Kornyshev (2017) Z. A. Goodwin and A. A. Kornyshev, Electrochemistry Communications 82, 129 (2017).
  * Goodwin  _et al._ (2017) Z. A. H. Goodwin, G. Feng, and A. A. Kornyshev, Electrochim. Acta 225, 190 (2017).
  * Ma  _et al._ (2015) K. Ma, J. Forsman, and C. E. Woodward, J. Chem. Phys. 142, 174704 (2015).
  * Lee  _et al._ (2014) A. A. Lee, D. Vella, S. Perkin, and A. Goriely, The journal of physical chemistry letters 6, 159 (2014).
  * Avni  _et al._ (2020) Y. Avni, R. M. Adar, and D. Andelman, Physical Review E 101, 010601 (2020).
  * Pivnic  _et al._ (2019) K. Pivnic, F. Bresme, A. A. Kornyshev, and M. Urbakh, Langmuir 35, 15410 (2019).
  * Choi  _et al._ (2017) J.-H. Choi, H. R. Choi, J. Jeon, and M. Cho, The Journal of chemical physics 147, 154107 (2017).
  * McEldrew  _et al._ (2020b) M. McEldrew, Z. A. Goodwin, S. Bi, M. Z. Bazant, and A. A. Kornyshev, J. Chem. Phys. 152, 234506 (2020b).
  * Tanaka (1989) F. Tanaka, Macromolecules 22, 1988 (1989).
  * Tanaka (1990) F. Tanaka, Macromolecules 23, 3784 (1990).
  * Tanaka and Stockmayer (1994) F. Tanaka and W. H. Stockmayer, Macromolecules 27, 3943 (1994).
  * Tanaka and Ishida (1995) F. Tanaka and M. Ishida, Journal of the Chemical Society, Faraday Transactions 91, 2663 (1995).
  * Ishida and Tanaka (1997) M. Ishida and F. Tanaka, Macromolecules 30, 3900 (1997).
  * Tanaka (1998) F. Tanaka, Physica A: Statistical Mechanics and its Applications 257, 245 (1998).
  * Tanaka and Ishida (1999) F. Tanaka and M. Ishida, Macromolecules 32, 1271 (1999).
  * Tanaka (2002) F. Tanaka, Polymer journal 34, 479 (2002).
  * Flory (1941) P. J. Flory, Journal of the American Chemical Society 63, 3083 (1941).
  * Stockmayer (1943) W. H. Stockmayer, The Journal of chemical physics 11, 45 (1943).
  * Eigen and Tamm (1962) M. Eigen and K. Tamm, Zeitschrift für Elektrochemie, Berichte der Bunsengesellschaft für Physikalische Chemie 66, 93 (1962).
  * Seo  _et al._ (2012a) D. M. Seo, O. Borodin, S.-D. Han, Q. Ly, P. D. Boyle, and W. A. Henderson, Journal of The Electrochemical Society 159, A553 (2012a).
  * Seo  _et al._ (2012b) D. M. Seo, O. Borodin, S.-D. Han, P. D. Boyle, and W. A. Henderson, Journal of The Electrochemical Society 159, A1489 (2012b).
  * Dutcher  _et al._ (2011) C. S. Dutcher, X. Ge, A. S. Wexler, and S. L. Clegg, The Journal of Physical Chemistry C 115, 16474 (2011).
  * Flory (1953) P. J. Flory,  _Principles of polymer chemistry_ (Cornell University Press, 1953).
  * Winter (2002) H. H. Winter, in  _Structure and Dynamics of Polymer and Colloidal Systems_ (Springer, 2002) pp. 439–470.
  * Stockmayer (1952) W. H. Stockmayer, Journal of Polymer Science 9, 69 (1952).
  * Christensen and Moloney (2005) K. Christensen and N. R. Moloney,  _Complexity and Criticality_ (PUBLISHED BY IMPERIAL COLLEGE PRESS AND DISTRIBUTED BY WORLD SCIENTIFIC PUBLISHING CO., 2005) https://www.worldscientific.com/doi/pdf/10.1142/p365 .
  * Choi  _et al._ (2018b) J.-H. Choi, H. Lee, H. R. Choi, and M. Cho, Annual review of physical chemistry 69, 125 (2018b).
  * Zheng  _et al._ (2018) J. Zheng, G. Tan, P. Shan, T. Liu, J. Hu, Y. Feng, L. Yang, M. Zhang, Z. Chen, Y. Lin, J. Lu, J. C. Neuefeind, Y. Ren, K. Amine, L.-W. Wang, K. Xu, and F. Pan, Chem 4, 2872 (2018).
  * Gonzales  _et al._ (2020) M. A. Gonzales, O. Borodin, M. Kofu, K. Shibata, T. Yamada, O. Yamamuro, K. Xu, D. L. Price, and M.-L. Saboungi, The Journal of Physical Chemistry Letters 11, 7279 (2020).
  * Garboczi  _et al._ (1995) E. Garboczi, K. Snyder, J. Douglas, and M. Thorpe, Physical review E 52, 819 (1995).
  * Ding  _et al._ (2017) M. S. Ding, A. von Cresce, and K. Xu, The Journal of Physical Chemistry C 121, 2149 (2017).
  * Ding and Xu (2018) M. S. Ding and K. Xu, The Journal of Physical Chemistry C 122, 16624 (2018).
  * Horwitz  _et al._ (2020) G. Horwitz, C. R. Rodríguez, P. Y. Steinberg, G. Burton, and H. R. Corti, Electrochimica Acta 359, 136915 (2020).
  * Steinrück  _et al._ (2020) H.-G. Steinrück, C. Cao, M. R. Lukatskaya, C. J. Takacs, G. Wan, D. G. Mackanic, Y. Tsao, J. Zhao, B. A. Helms, K. Xu, O. Borodin, J. F. Wishart, and M. F. Toney, Angewandte Chemie (2020).
  * Borodin  _et al._ (2017b) O. Borodin, X. Ren, J. Vatamanu, A. von Wald Cresce, J. Knap, and K. Xu, Accounts of chemical research 50, 2886 (2017b).
  * Steinrück  _et al._ (2018) H.-G. Steinrück, C. Cao, Y. Tsao, C. J. Takacs, O. Konovalov, J. Vatamanu, O. Borodin, and M. F. Toney, Energy & Environmental Science 11, 594 (2018).
  * Kohns  _et al._ (2016) M. Kohns, S. Reiser, M. Horsch, and H. Hasse, The Journal of chemical physics 144, 084112 (2016).
  * Xu  _et al._ (2007) K. Xu, Y. Lam, S. S. Zhang, T. R. Jow, and T. B. Curtis, The Journal of Physical Chemistry C 111, 7411 (2007).
  * Xu (2007) K. Xu, Journal of The Electrochemical Society 154, A162 (2007).
  * von Wald Cresce  _et al._ (2012) A. von Wald Cresce, O. Borodin, and K. Xu, The Journal of Physical Chemistry C 116, 26111 (2012).
  * Zhang  _et al._ (2018) H. Zhang, B. Qin, J. Han, and S. Passerini, ACS Energy Letters 3, 1769 (2018).
  * Dou  _et al._ (2018) 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. 11, 3212 (2018).
  * Molinari and Kozinsky (2020) N. Molinari and B. Kozinsky, J. Phys. Chem. B 124, 2676 (2020).
  * Lui  _et al._ (2011) M. Y. Lui, L. Crowhurst, J. P. Hallett, P. A. Hunt, H. Niedermeyer, and T. Welton, Chem. Sci. 2, 1491 (2011).
  * Suo  _et al._ (2016) L. Suo, O. Borodin, W. Sun, X. Fan, C. Yang, F. Wang, T. Gao, Z. Ma, M. Schroeder, A. von Cresce, S. M. Russell, M. Armand, A. Angell, K. Xu, and C. Wang, Angewandte Chemie International Edition 55, 7136 (2016).
  * Kondou  _et al._ (2018) S. Kondou, E. Nozaki, S. Terada, M. L. Thomas, K. Ueno, Y. Umebayashi, K. Dokko, and M. Watanabe, The Journal of Physical Chemistry C 122, 20167 (2018).
  * Molinari  _et al._ (2019a) N. Molinari, J. P. Mailoa, N. Craig, J. Christensen, and B. Kozinsky, Journal of Power Sources 428, 27 (2019a).
  * Molinari  _et al._ (2019b) N. Molinari, J. P. Mailoa, and B. Kozinsky, The journal of physical chemistry letters 10, 2313 (2019b).
  * Chen  _et al._ (2020b) L. Chen, J. Zhang, Q. Li, J. Vatamanu, X. Ji, T. P. Pollard, C. Cui, S. Hou, J. Chen, C. Yang, L. Ma, M. S. Ding, M. Garaga, S. Greenbaum, H.-S. Lee, O. Borodin, K. Xu, and C. Wang, ACS Energy Letters (2020b).
  * Jiang  _et al._ (2020) L. Jiang, L. Liu, J. Yue, Q. Zhang, A. Zhou, O. Borodin, L. Suo, H. Li, L. Chen, K. Xu, and Y. Hu, ADV MATER 32, 1904427 (2020).
  * Becker  _et al._ (2020) M. Becker, D. Reber, A. Aribia, C. Battaglia, and R.-S. Kühnel, (2020).
  * Flory (1942) P. J. Flory, The Journal of chemical physics 10, 51 (1942).
  * Flory (1956) P.-J. Flory, Proceedings of the Royal Society of London. Series A. Mathematical and Physical Sciences 234, 60 (1956).
  * Tanaka (2011) F. Tanaka,  _Polymer physics: applications to molecular association and thermoreversible gelation_ (Cambridge University Press, 2011).
  * Plimpton (1995) S. Plimpton, J. Comput. Phys. 117, 1 (1995).
  * Martínez  _et al._ (2009) L. Martínez, R. Andrade, E. G. Birgin, and J. M. Martínez, J. Comput. Chem. 30, 2157 (2009).
  * Humphrey  _et al._ (1996) W. Humphrey, A. Dalke, and K. Schulten, J. Mol. Graphics 14, 33 (1996).
  * Lopes and Pádua (2012) J. N. C. Lopes and A. A. Pádua, Theoretical Chemistry Accounts 131, 1129 (2012).