# Cluster Formation induced by local dielectric saturation in Restricted
Primitive Model Electrolytes

David Ribar  Computational Chemistry, Lund University, P.O.Box 124, S-221 00
Lund, Sweden  Clifford E. Woodward  School of Physical, Environmental and
Mathematical Sciences University College, University of New South Wales, ADFA
Canberra ACT 2600, Australia  Sture Nordholm  Department of Chemistry and
Molecular Biology, The University of Gothenburg, 412 96 Gothenburg, Sweden
Jan Forsman<EMAIL_ADDRESS>Computational Chemistry, Lund
University, P.O.Box 124, S-221 00 Lund, Sweden

###### Abstract

Experiments using the Surface Force Apparatus (SFA) have found anomalously
long ranged charge-charge underscreening in concentrated salt solutions.
Meanwhile, theory and simulations have suggested ion clustering to be the
possible origin of this behaviour. The popular Restricted Primitive Model of
electrolyte solutions, in which the solvent is represented by a uniform
relative dielectric constant, $\varepsilon_{r}$, is unable to resolve the
anomalous underscreening seen in experiments. In this work, we modify the
Restricted Primitive Model to account for local dielectric saturation within
the ion hydration shell. The dielectric constant in our model locally
decreases from the bulk value to a lower saturated value at the ionic surface.
The parameters for the model are deduced so that typical salt solubilities are
obtained. Our simulations for both bulk and slit geometries show that our
model displays strong cluster formation and these give rise to long-ranged
interactions between charged surfaces at distances similar to what has been
observed in SFA measurements. An electrolyte model wherein the dielectric
constant remains uniform does not display similar clusters, even with
$\varepsilon_{r}$ equal to the saturated value at ion contact. Hence, the
observed behaviours are not simply due to an enhanced Coulomb interaction.

###### keywords:

Restricted Primitive Model, Electrolytes, Dielectric saturation, Simulations,
Anomalous screening, Ion clustering

## 1

Electrolytes play an important role in a plethora of both scientific and
industrial applications 1, 2, 3. Simple theoretical descriptions at the mean-
field level 4, 5, 1 have proven to be reasonably accurate for aqueous
electrolytes at low coupling strength, e.g., monovalent ions (1:1 salts) in
aqueous solvents at low and intermediate concentrations. Perhaps the most
fundamental prediction by these theories is so-called ionic screening, often
described in terms of the so-called Debye screening length, $\lambda_{D}$. The
Debye length is predicted to be inversely proportional to the square root of
the electrolyte concentration, $\lambda_{D}\sim 1/c^{1/2}$. Salts composed of
ions where at least one component is multivalent will generally require a
higher level of theory, accounting for ion-ion correlations 6, 7, 8, 9. This
notwithstanding, it was believed that mean-field theories could capture the
qualitative behaviours of even concentrated aqueous solutions of 1:1 salts,
with some predictable corrections due to correlations. This has been called
into question by recent experimental investigations of 1:1 electrolyte
solutions, above a threshold concentration of about 1 M. These have revealed a
peculiar anomalous underscreening phenomenon 10, 11, 12, 13, 14, 15 whereby
the interaction between charged surfaces exhibit an exponential decay with an
extraordinarily long range, $\lambda$, much larger than the Debye length
predicted by mean-field theory. Moreover, $\lambda$ appears to increase with
the salt concentration, in qualitative disagreement with more sophisticated
theories that attempt to correct for correlations 16, 17, 18. It should be
noted that anomalous underscreening has been experimentally challenged 19, and
despite considerable theoretical efforts 17, 20, 21, 22, 23, 24, 18 there is
at present no consensus as to its physical origin.

One often proposed explanation involves the role of clusters and their effect
on the Debye screening length. For example, in a concentrated 1:1 electrolyte,
the formation of neutral and/or weakly charged clusters in concentrated ionic
solutions would act to reduce the number of independent charged species. If,
$c*$, denotes the concentration of the effective screening charge (and we
assume most clusters are either neutral or univalent), then the modified Debye
length would be given by, $\lambda\sim 1/c*^{1/2}$. Screening length
measurements in the anomalous region have been fit to the relation,
$\lambda\sim cl_{B}d^{3}$, where $l_{B}=\beta
e_{0}^{2}/(4\pi\varepsilon_{0}\varepsilon_{r})$ is the Bjerrum length and $d$
is an average ionic diameter. Here $\beta=(k_{B}T)^{-1}$ is the inverse
thermal temperature, $e_{0}$ denotes the elementary charge, and
$\varepsilon_{0}$ is the the permittivity of vacuum. Attributing
underscreening to clustering and a modified Debye length, would require
$c^{*}\sim 1/c^{2}$ in the anomalous region.

It is worthwhile to reflect upon the implications of such an explanation in a
little more detail. Suppose, we consider an electrolyte solution with a
concentration well below the threshold value ($\sim$ 1 M) above which
underscreening occurs. Presumably, there will be some incipient clustering of
ions at this concentration, but insufficient to affect the screening length
significantly. Clustering at low concentration is driven by the electrostatic
attractions between ions, but clusters remain finite due to the low chemical
potential of the free ions (with which the aggregated ions remain in
equilibrium). As the concentration increases, clustering is additionally aided
by the lowering of the overall excluded volume due to ion aggregation. At
concentrations lower than the threshold value for underscreening, addition of
more ions will tend to lower the screening length, as the number of free ions
at equilibrium will generally increase with the overall concentration. Once
the concentration reaches the underscreening threshold value, however,
addition of more ions will cause the screening length to increase. Since we
are now in the region of anomalous underscreening, effectively all the added
ions will be aggregated to create neutral clusters. In addition, a portion of
the original population of charged species (charged clusters and free ions)
must also reorganise into neutral structures. If correct, this has the
hallmarks of an apparent instability in the free energy of cluster formation.
It suggests that, at a threshold concentration of free ions, clusters become
unstable and experience accelerated growth, presumably because they have
surpassed a critical size. This growth is arrested by a concurrent decrease in
the free ion concentration. Experiments indicate that the concentration of
effective screening charges would need to be some $10^{4}$ times lower than
$c$ in order to explain the upper levels of underscreening 10, 25, 11, 12.

In electrolyte theory the Restricted Primitive Model (RPM) has been something
of a "workhorse". Here ions are treated as charged hard spheres with diameter,
$d$, and the solvent is simply modelled using a uniform dielectric constant,
$\varepsilon_{r}$, which is generally greater than 1 due to the polarisability
of solvent molecules. Thus, ion-ion interactions are described by,
$\phi_{ij}(r)$, where:

$\beta\phi_{ij}(r)=\left\\{\begin{array}[]{ll}\infty;&r\leq d\\\
l_{B}\frac{z_{i}z_{j}}{r};&r>d\\\ \end{array}\right.$ (1)

Here, $z_{i}$ and $z_{j}$ denote the valencies of interacting ions $i$ and
$j$. Interestingly, hard sphere interactions tend to shorten the electrostatic
screening length compared with $\lambda_{D}$, due to interplay between hard
sphere and electrostatic correlations. 16.

A clustering mechanism was recently explored using simulations of the RPM18.
The simulations showed that the propensity of clusters to form increases when
the coupling strength and/or the electrolyte concentration increases and that
this does lead to underscreening. Furthermore, a clustering analysis was used
to confirm that the measured screening length was consistent with using the
apparent concentration of independent charge carriers, $c*$ (charged clusters
and free ions), i.e., $\lambda\sim 1/c*^{1/2}$. However, the underscreening
observed was insufficient to explain experimental results. In particular, the
simulations did not predict an increase in the observed screening length as
the overall electrolyte concentration increases, i.e., the simulated $c*$
always increased with $c$. It should be noted that these simulations also
accounted for a reduction in $\varepsilon_{r}$ with increasing $c$. Thus the
RPM is not able to reproduce the trends seen in experiments, which is perhaps
not unanticipated given the discussion above. That is, the RPM appears to have
no inherent features that would suggest the possibility of anomalous
clustering of the type discussed above.

Given this, it is interesting to consider other explanations (apart from
clustering) to explain anomalous underscreening. In particular, we consider
here the possibility that the decay length in experiments may be a reflection
of the size of clusters, rather than their impact on the concentration of free
charges. In the case of NaCl, the largest decay lengths measured in the
anomalous regime are of the order of $30$ Å. 12 A pair of $Na^{+}$ $Cl^{-}$
ions in contact has a size of $\sim 6$ Å, while water molecules have a
diameter $\sim 3$ Å. Thus, one could envisage a small cluster of a few pairs
of ions and hydrating water molecules to have a size of at least $\sim 30$ Å.
While such a mechanism has not been uncovered for the RPM, that may be a
consequence of short-comings in that interaction model, which does not allow
the formation of large enough clusters.

Most simple theories and simulation models of electrolytes do not explicitly
include the solvent but allows it to be represented in a “primitive” fashion,
via a relative dielectric constant, $\varepsilon_{r}$, which is greater than
unity due to the interaction of the ion charge with electrons and nuclear
charges of the surrounding solvent molecules. Water is a polar solvent and can
reorient and rearrange in response to the electric field of an ion.
Experimentally, it is known that as the concentration of an aqueous
electrolyte solution increases the overall dielectric constant decreases 26,
27, 28, 29, 30, 31. This can be rationalised if we consider that hydrating
water molecules are rotationally constrained by electric fields from the ions.
As the salt concentration increases to molar levels, these constrained water
molecules constitute an increasingly significant percentage of the solvent,
which leads to a decreased overall dielectric response. This type of
dielectric saturation has been the focus of some previous investigations at
electrode surfaces, or in confined geometries29, 27, 32. Theoretical
investigations have explored this reduction in dielectric permittivity by
correcting the mean field approach in a consistent manner, 21 or else model it
in terms of ion-specific effects 33, 34. In protein small angle x-ray
scattering model fitting procedures, the formation of water hydration shells
around charged (and uncharged) proteins, is a well-established phenomenon. A
special treatment of hydration shells is needed, in order to produce a model
fit that correctly describes the measured spectra 35, 36, 37.

We propose to modify the RPM, based on physically plausible arguments. To
maintain simplicity we will explore a 2-body interaction model, but ultimately
the effects we consider are best manifested using many-body forces. This
notwithstanding, we consider it useful to understand the effect on screening
lengths in the presence of much larger clusters than can be generated by the
simple RPM. We use this model to roughly mimic NaCl and show that it saturates
at a concentration close to its observed value, if not somewhat lower. While
we by no means present this new model as an accurate representation of aqueous
NaCl solutions, it does allow us to investigate screening behaviour in
electrolytes in the presence of much larger clusters than is generated by the
RPM. This allows us to gain some further insight as to potential mechanisms
that a cluster model can provide to explain anomalous underscreening apart
from the usual assertion that it reduces the concentration of free charges.

A physically reasonable modification of the RPM, to account for dielectric
saturation, would be to make $\varepsilon_{r}$ in Eq.(1) a function of the
average electrolyte concentration, as has been employed in previous RPM
studies. 18. This will increase the coupling strength between ions at higher
average electrolyte concentration. However, even in dilute solutions, we
expect that dielectric saturation will occur in regions where fluctuations
cause a locally high concentration of ions. In such a region, solvent
molecules will be displaced, and those that remain are subject to large
electric fields 29. This suggests a spatial inhomogeneity in $\varepsilon_{r}$
should occur as a result of changes in the ionic configurations. In principle,
dielectric inhomogeneity should be modelled as a many-body phenomenon that
qualitatively results in local changes to the electrostatic interactions. A
collection of ions will produce large fields reducing the local dielectric
constant and promoting clustering. Countering this will be an energetically
unfavourable contribution due to the overlap of solvation shells of those
ions. In this work, we consider a phenomenological manifestation of these
mechanisms at the pair-wise interaction level, in order to maintain
computational simplicity.

In the model used here, $\varepsilon_{r}$, is assumed to vary with the ionic
separation, $r$. Specifically, it will have a reduced value when ions approach
each other but has its bulk value when the ions are sufficiently separated. We
will model $\varepsilon_{r}(r)$ as a linear ramp function of the following
form:

$\varepsilon_{r}(r)=\left\\{\begin{array}[]{ll}\varepsilon_{c};&r\leq d\\\
\varepsilon_{c}+(\varepsilon_{b}-\varepsilon_{c})\frac{r-d}{\Delta};&d<r\leq
d+\Delta\\\ \varepsilon_{b};&r>d+\Delta\\\ \end{array}\right.$ (2)

with $\varepsilon_{c},\varepsilon_{b}$ denoting the contact and bulk solvent
dielectric constant values, respectively. In general we have $\varepsilon_{c}$
< $\varepsilon_{b}$ The slope of the linear ramp is defined by the parameter
$\Delta$, which is taken as the diameter of a solvent molecule, $d$. We chose
a value of $d=$ 3 Å representing a single hydration layer of water molecules.
Consistent with this, we set $\varepsilon_{b}=78.3$. The thus modified
potential used here, consists essentially of a short-ranged interaction
(attractive between unlike ions, repulsive between like ions) in addition to a
typical RPM interaction, wherein the latter assumes the uniform dielectric
constant of the bulk solvent, $\varepsilon_{b}$. While other implementations
of the RPM have assumed a uniform dielectric constant, $\varepsilon_{r}$, that
decreases with electrolyte concentration to account dielectric saturation, we
shall see below that a locally dependent $\varepsilon_{r}$ of the type
proposed here promotes more clustering, even compared to an RPM that uses
$\varepsilon_{r}=\varepsilon_{c}$ everywhere. We will use our potential model
to explore both bulk simulations of the electrolyte as well as the solution in
contact with charged surfaces.

We will only give a brief account of the simulation methods here, and refer to
the Supplementary Information (SI) for details. Canonical ensemble Metropolis
Monte-Carlo simulations were performed at 298 K for two geometries: a cubic
simulation geometry, henceforth referred to as the bulk, and a parallelepiped
geometry between two charged hard walls, referred to as the slit. Periodic
boundary conditions were applied along $(x,y,z)$ for the bulk simulations
(with side-length of the cubic simulation box denoted as $L$) and along
$(x,y)$ for the slit simulations, with impenetrable charged hard walls
situated at $z=\pm H/2$. For slit systems, $L$ denotes the length of the
simulation parallelepiped along the $(x,y)$-axis, and $H$ represents the slit
width along the $z$-axis. The values of $L$ and $H$ are adjusted to
accommodate the different electrolyte concentrations under investigation. For
the bulk simulations, we employed Minimum Image (MI) truncation of the Coulomb
interactions. The accuracy of this choice has been evaluated thoroughly in a
recent work 38, by direct comparisons with more elaborate (and computationally
expensive) Ewald simulations. These comparisons verified that MI truncation
leads to structurally accurate results for bulk systems. We have included yet
another comparison, with the same conclusion, in the SI. For the slit systems,
we employed the standard “charged sheet” method 39 to account for long-ranged
interactions. Notably, cluster moves were implemented for both systems,
leading to crucial improvements of the statistical performance.

There are two free parameters in our model potential, which control the
magnitude of the short-ranged attraction. We found that at a high coupling
strength, (small $d$ and $\varepsilon_{c}$) the bulk solution appears to
separate into an amorphous condensed phase in equilibrium with a dilute
clustered phase. This could be viewed as the saturation limit for the
solution, except in real electrolytes the condensed phase would be an ordered
crystal. This phase instability is seen in Figure 1(a) where we observe the
response of the cation-cation pair correlation $g_{++}(r)$ to changes of
$\varepsilon_{c}$, at a concentration of 3.45 M. The development of a
pronounced long-ranged slope in the tail of $g_{++}(r)$ is a signature of a
phase separation. This was further supported by configurational snapshots in
the SI, whereby visual inspection of the condensed phase suggests that it is
non-crystalline in nature. While a more thorough structural analysis is
lacking, we assert that determining the precise nature of the condensed phase
is anyway not pertinent to the aims of this study. We can, however, plausibly
argue that an amorphous condensed phase is not surprising, given the
approximate nature of our model, wherein significant many-body effects are
ignored. In particular, many-body effects are expected to promote clusters
with ordered cores but labile outer regions, as dielectric saturation is
larger towards the centre of the clusters, while closer to the extremities
solvent polarisation remains large. On the other hand, the pair-wise
approximation, used in this study, somewhat erroneously strengthens ion-ion
interactions throughout the cluster, as ameliorated by the choice of
$\varepsilon_{c}$. Thus, in order to predict a reasonable saturation
concentration for the electrolyte model, a sensible choice for
$\varepsilon_{c}$ should fall between the bulk value of the solvent and the
(small) internal dielectric constant of a crystal. However, such a choice will
likely stabilise an amorphous condensed phase rather than a crystalline phase,
given that $\varepsilon_{c}$ will be larger than the expected dielectric
response inside a crystal. Furthermore, the amorphous phase will also have
less binding energy than that of the crystal, which suggests that our pair-
wise approximation will possibly predict less clustering as the solution
concentration approaches saturation. KCl and NaCl have saturation
concentrations of about 4.5-6 M, which suggests we should consider a similar
saturation concentration in our modelling.

Given that the aim of our study was primarily to explore the impact of
clustering on electrostatic correlations, it was important that we select a
value for $\varepsilon_{c}$ so that the system displayed large but finite
clusters. From Figures 1(a)-(b) we see that, at 3.45 M, the system phase
separates for $\varepsilon_{c}$ values of 20 and 21, but not when
$\varepsilon_{c}$ is above 22 (for $d=$ 3 Å) . Thus we have set
$\varepsilon_{c}=23$ for all subsequent simulations. This value provides a
significant degree of clustering in our model while the solution remains
unsaturated at least up to 3.45 M. A different choice of $d$ would lead to a
different choice of $\varepsilon_{c}$, as described in the Supporting
Information, SI

(a)

(b)

Figure 1: (a) Effect of changing the $\varepsilon_{c}$ on the phase stability
at 3.45 M. All systems contain 5000 ion pairs. (b) An illustration of the
progression from a phase separated to homogenised system, for a system with
$d=3$ Å and $\varepsilon_{c}=23$. The phase separated system was created from
simulations with $\varepsilon_{c}=21$. Displayed are the cation-cation radial
distribution functions from a set of short simulations (dashed lines), along
with a sufficiently long final equilibrium simulation (with
$\varepsilon_{c}=23$).

It is of course possible that the system is metastable rather than stable for
$\varepsilon_{c}=23$ at 3.45 M. However, we have made tests ensuring that a
metastable scenario is highly unlikely. Specifically, we have (several times)
initiated simulations from a phase separated system, at $\varepsilon_{c}=21$,
only to find that, upon switching to $\varepsilon_{c}=23$, the system
transitions to a single phase state at equilibrium. This is illustrated in
Figure 1(b), where we depict how the pronounced long-ranged gradient of the
initial radial distribution function (indicating a phase-separated system),
gradually disappears.

Figure 2: Results from bulk simulations. (a) Cation-cation (or anion-anion)
resolved radial distribution functions, $g_{++}(r)$ (on average identical to
$g_{--}(r)$), for the local dielectric saturation model (full line), and a
uniform low-dielectric constant (a uniform value of $\varepsilon_{r}=23$)
model (dashed line). (In order to improve statistics we have in reality
measured “$g_{++}(r)$” as $(g_{++}(r)+g_{--}(r))/2$) (b) Cation-anion resolved
radial distribution functions, $g_{+-}(r)$. (c) $g_{cc}(r)\equiv
g_{+-}(r)-(g_{++}(r)+g_{--}(r))/2$. (d) Probability of an ion to be part of a
cluster of size $N_{c}$, $P_{ion}(N_{c})$, at three different concentrations.
An ion must be within a distance $\delta$ or less, of at least one other ion
within a cluster, in order to be a member of that cluster. We have set
$\delta=d+0.5$Å. The inset is a focus on small clusters.

We reiterate that our potential model introduces an additional short-ranged
asymmetric interaction, which promotes association between unlike ions and the
formation of large clusters. This cluster formation as a function of the
electrolyte concentration can be readily observed in the species resolved
radial distribution functions, $g_{ij}(r)$ (with $i,j=\pm$) as illustrated in
Figure 2(a)-(b). By symmetry $g_{--}(r)=g_{++}(r)$ and $g_{+-}(r)=g_{-+}(r)$
thus it is prudent to use $(g_{++}(r)+g_{--}(r))/2$ and
$(g_{+-}(r)+g_{-+}(r))/2$ for the like and unlike correlation functions. We
can define the total density correlation function
$g_{nn}(r)=(g_{++}(r)+g_{--}(r)+g_{+-}(r)+g_{-+}(r))/4$ (the particle density
correlation around any ion) as well as the so-called charge-charge correlation
function $g_{cc}(r)=(g_{+-}(r)+g_{-+}(r)-g_{++}(r)-g_{--}(r))/2$ (the counter-
charge density around an ion). The charge-charge correlation functions are
shown in Figure 2(c). Using linear response theory it is possible to relate
the interaction free energy between two charged surfaces (at large separation)
to these correlation functions at long-range. As charged surfaces perturb both
the charge and particle densities of the electrolyte contained between them,
the interaction between the surfaces at large separations is dictated either
by $g_{cc}(r)$ or $g_{nn}(r)$, whichever has the longest range. Classical
mean-field theory, wherein ions are assumed to respond only to the mean
electrostatic potential, asserts that $h_{cc}(r)$ has a Yukawa form,
$h_{cc}(r)\sim\exp(-r/\lambda_{D})/r$ and $h_{nn}(r)$ is much shorter-ranged
($h_{\alpha\beta}(r)=g(r)_{\alpha\beta}-1$). In fact, linearised mean-field
analysis predicts that the charged surfaces will not alter the total ionic
density between them at all, $h_{nn}(r)=0$, whereas to second order we have
$h_{nn}(r)\sim\exp(-2r/\lambda_{D})/r^{2}$. This qualitative picture may
change in the presence of strong correlations and, as in our case, with the
formation of large clusters.

At the lowest concentration investigated (0.05 M), we observe what could be
described as mean-field behaviour for the correlation functions. Here
$g_{++}(r)$ displays a co-ion exclusion region, while $g_{+-}(r)$ illustrates
counter-ion attraction, so that $h_{++}(r)\approx-h_{+-}(r)$ as they approach
zero. Between 0.25-0.5 M a peak appears in $g_{++}(r)$ at short-range, which
then displays co-ion exclusion at larger distances. As the concentration
increases however, the co-ion exclusion region diminishes in size so that at
0.5 M it has almost vanished. On the other hand, $g_{+-}(r)$ still displays an
attraction between counterions, but its range appears to decrease, which is
actually qualitatively consistent with the mean-field behaviour. In Figure
2(c), we observe the a reduction in range of $g_{cc}(r)$ as well. In addition,
oscillations appear at short range at 0.5 M, which, together with the
behaviour of $g_{++}(r)$, suggests accumulation of ions into clusters.
However, even in this concentration range, the tails of these correlation
functions are still Yukawa like and as we show below (and in the SI)
asymptotic fitting gives a correlation length which is greater than
$\lambda_{D}$, indicating some degree of underscreening.

At even higher concentrations (1 M and 3 M), we note the complete
disappearance of the co-ion exclusion region in $g_{++}(r)$. This is
accompanied by apparent coincidence of $g_{++}(r)$ and $g_{+-}(r)$ at longer
range, as indicated by the diminished range of their difference, $g_{cc}(r)$,
as shown in Figure 2(c). That is, at lower concentrations, the long-range
decays of all these correlation functions were similar, albeit $g_{++}(r)$ and
$g_{+-}(r)$ approach unity from below and above respectively. Indeed, in the
SI we show that they can be all reasonably well fitted to a Yukawa form,
$\sim\exp(-r/\lambda)/r$, where $\lambda$ is generally larger than the Debye
length. However, at 1 M and 3 M, $g_{++}(r)$ and $g_{+-}(r)$ both approach
unity from above and, their range is significantly longer than $g_{cc}(r)$.
These results indicate that there is a great deal of clustering occurring,
causing much longer range correlations in the density, $g_{nn}(r)$, compared
with the charge, $g_{cc}(r)$. The correlation length of $g_{nn}(r)$ will be
determined by the typical cluster size. In Figure 2(d), we show the
probability, $P_{ion}(N_{c})$, that an ion is a member of a cluster of size
$N_{c}$ for different concentrations, and the growth in cluster sizes is
apparent as the ion concentration increases. The growth in clusters is
accompanied with more rapid screening of charges. This effect on charge
screening is also qualitatively predicted by mean-field theories, as the Debye
length decreases with concentration. But here we also see that $g_{cc}(r)$
loses its Yukawa form and becomes oscillatory, so hard-core correlations
within dense clusters are clearly playing a role. Most clusters are expected
to be close to neutral, which is why the difference between cation-cation and
cation-anion correlations at long range vanishes. What has not been previously
reported in other theoretical treatments of electrolyte models (as far as we
are aware), is the significant increase in range of $g_{nn}(r)$ as a function
of concentration, compared with $g_{cc}(r)$.

Thus, our model predicts that at around 1 M concentration, the interaction
between charged surfaces will begin to be dominated by density correlations,
as measured by $g_{nn}(r)$, rather than charge correlations (from
$g_{cc}(r)$), which were significant at lower concentrations. In our model,
density correlations are affected by cluster formation and their typical size,
whereas charge correlations are determined by the availability of screening
charges.

We argue that our modification of the RPM gives rise to clusters that are much
larger than usual applications of the RPM. To illustrate this, we also
calculated the correlation functions at 3 M for a simple RPM case, where the
uniform dielectric constant was chosen to be
$\varepsilon_{r}=\varepsilon_{c}=23$. That is, we assumed dielectric
saturation occurs over the full range of the Coulomb interaction, rather than
just close to ion-ion contact. We see (Figure 2(a)) that $g_{++}(r)$ then
displays a short-ranged co-ion exclusion region superimposed on an oscillatory
profile, with no long-ranged tail as observed in our modified RPM. Similarly,
$g_{+-}(r)$ shows a short-ranged counter-ion enhancement with oscillations,
again without a long-ranged tail (Figure 2(b) . The corresponding $g_{cc}(r)$
is also oscillatory, but with amplitude much smaller than the modified RPM
(Figure 2(c)), which is evidence of much larger clustering induced by the
modified RPM. It is clear that the additional short-ranged interaction
introduced in our modified RPM reduces the incentive for ions to dissociate
from clusters once they are formed. If dielectric saturation is assumed to
occur everywhere, even in regions of low ionic density, there is less of an
advantage for ions to cluster. Note, that this affect is ameliorated by the
decreased in ion solvation expected to occur in clusters, an effect which must
be accounted for when choosing a suitable value for $\varepsilon_{c}$.

It may seem surprising that the range of a reduced dielectric response has
such a strong impact on the cluster forming tendency. However, one should take
cognisance of the fact that with a uniform and low dielectric constant, there
are inevitable strong repulsions between like charges in a cluster. These
repulsions are considerably weaker if the reduced dielectric response is
local, mainly influencing the interaction between charges of opposite sign.

Figure 3: Effective electrostatic screening lengths obtained via the modified
Widom technique, and the corresponding Debye screening lengths. The dashed
lines guide the eye.

In previous work, we developed a simulation method that uses a modified Widom
technique to estimate the long range decay length of $g_{cc}(r)$, which is
very accurate if one can assume a Yukawa form, $g_{cc}(r)\sim
1+A\exp(-\kappa_{eff}r)/r$ 38. The method is briefly described in the SI. The
results for the correlation functions calculated above are compared with the
standard Debye screening length ($\kappa_{D}^{2}=\beta c\sum_{i}(ez_{i})^{2}$)
in Figure 3. Note that while the Yukawa form is only really valid below 1 M,
our method is nevertheless still able to extract an effective screening length
at the higher concentrations. This is because it is based on a linear
expansion of the free energy functional that predicts a screening length from
an ensemble average. This average is still calculable even in cases where
assumptions which lead to a Yukawa form break down. In Figure 3 we do observe
some underscreening, the degree of which is much smaller than that suggested
by SFA experiments10, 11, 12, 13, 14. In particular, we do not observe an
increase in screening length with concentration. In any case, as described
above, the decay length of $g_{cc}(r)$ is not the dominant one above 1 M in
any case, but rather that of $g_{nn}(r)$. Interestingly, by the use of
asymptotic analysis techniques, researchers have previously determined that
such a scenario is theoretically possible 40, 16. That is, even in the RPM, a
transition may occur from charge-charge correlation to a density-density
dominated correlation in the asymptotic (long-range) regime.

Figure 4: Results from the slit geometry simulations, with negatively charged
surfaces. (a) Concentrations of cations (solid lines) and anions (dashed
lines) along the $z$-axis, for three different simulated bulk concentrations
(these are estimated from mid plane values). (b) The concentration difference
between cations and anions along the $z$-axis. Here, we have also added data
from RPM simulations with a uniform dielectric constant of 23, at 1.4M (dashed
line).

Finally, we turn our attention to structural properties of our modified RPM,
in the presence of two macroscopic flat and negatively charged surfaces. For
these systems, we have chosen a simulation model with a slit geometry. The
slit simulations enable us to study the behaviour of the model in a system
that approximates the experimental SFA setup: an electrolyte confined between
two macroscopic charged surfaces. We investigated three concentrations at a
constant inverse surface charge density of $-70\,e/$Å2. The resulting ion
density distributions, $n_{+}(z)$ and $n_{-}(z)$, where $z$ is the direction
normal to the surfaces, are presented in Figure 4(a). We can observe a
significant difference between $n_{+}(z)$ and $n_{-}(z)$ at the mid-plane
between the surfaces, even when these are 50 Å apart. This is in stark
contrast to mean-field predictions, since the mid plane is more than 13 Debye
lengths distant from the surfaces, at the highest investigated bulk
concentration (about 2.6 M). In order to clarify our results this further, we
plot the density difference, $\Delta n(z)\equiv n_{+}(z)-n_{-}(z)$, in Figure
4(b). Even though the overall $\Delta n(z)$ profile does drop as the salt
concentration increases, we note that the salt dependence is quite weak, and
that a significant long-ranged tail persists also for very high ionic
strengths. The scenario would again be quite different with a model using a
uniform dielectric constant (even $\varepsilon_{r}$=23), in which case $\Delta
n(z)$ rapidly vanishes away from a charged surface, at high concentrations.
This is explicitly shown by the dashed grey line in Figure 4(b). While we have
not measured the net interaction between such charged surfaces in this study,
the long-ranged tail of $\Delta n(z)$ clearly implies a slowly decaying
surface force. This force will be quantified in future simulation work.

Considering the fact that experimental approaches extract the effective
correlation lengths14, which may or may not be purely a result of charge-
charge correlations, the observation of a transition from the expected charge-
charge asymptotic domination to a density-density domination at high
concentrations strongly supports the hypothesis that cluster formation, has a
crucial influence on anomalous underscreening effects.

This theoretical study has investigated the influence of local dielectric
saturation on the structure of electrolytes modelled by the RPM. We
demonstrate that local dielectric saturation induces significant clustering.
Interactions between such clusters dominate the asymptotic correlations for
systems approaching or exceeding concentrations of about 1 M. When charged
surfaces are immersed in such solutions, a net charge density develops, that
decays quite slowly with the transverse distance to these surfaces. This
suggests that there might be long-ranged interactions between such surfaces,
commensurate with observations by SFA.

J.F. acknowledges financial support by the Swedish Research Council, and
computational resources by the Lund University computer cluster organisation,
LUNARC.

The following files are available free of charge.

  * •

Supporting information: Detailed simulation methods, and further analyses.

  * •

Github repository: all codes used for simulations, along with the data
generated, is freely available.

## References

  * Israelachvili 1991 Israelachvili, J. N. _Intermolecular and Surface Forces, 2nd Ed._ ; Academic Press: London, 1991
  * Evans and Wennerström 1994 Evans, F. A.; Wennerström, H. _The colloidal domain: where Physics, Chemistry, Biology and Technology meet_ ; VCH Publishers: New York, 1994
  * Holm et al. 2001 Holm, C.; Kekicheff, P.; Podgornik, R. _Electrostatic Effects in Soft Matter and Biophysics_ ; Kluwer Academic Publishers: Dordrecht, 2001
  * Derjaguin and Landau 1941 Derjaguin, B. V.; Landau, L. Theory of the Stability of Strongly Charged Lyophobic Sols and of the Adhesion of Strongly Charged Particles in Solutions of Electrolytes. _Acta Phys. Chim. URSS_ 1941, _14_ , 633–662
  * Verwey and Overbeek 1948 Verwey, E. J. W.; Overbeek, J. T. G. _Theory of the Stability of Lyophobic Colloids_ ; Elsevier Publishing Company Inc.: Amsterdam, 1948
  * Nordholm 1984 Nordholm, S. Generalized van der Waals theory. XII. Application to ionic solutions. _Aust. J. Chem._ 1984, _37_ , 1
  * Guldbrand et al. 1984 Guldbrand, L.; Jönsson, B.; Wennerström, H.; Linse, P. Electrical Double Layer Forces, A Monte Carlo Study. _J. Chem. Phys._ 1984, _80_ , 2221
  * Kjellander and Marcelja 1986 Kjellander, R.; Marcelja, S. Interaction of charged surfaces in electrolyte solutions. _Chem. Phys. Lett._ 1986, _127_ , 402–407
  * Valleau et al. 1991 Valleau, J.; Ivkov, R.; Torrie, G. M. _J. Phys. Chem._ 1991, _95_ , 520
  * Gebbie et al. 2013 Gebbie, M. A.; Valtiner, M.; Banquy, X.; Fox, E. T.; Henderson, W. A.; Israelachvili, J. N. Ionic liquids behave as dilute electrolyte solutions. _PNAS_ 2013, _110_ , 9674–9679
  * Gebbie et al. 2015 Gebbie, M. A.; Dobbs, H. A.; Valtiner, M.; Israelachvili, J. N. Long-range electrostatic screening in ionic liquids. _PNAS_ 2015, _112_ , 7432–7437
  * Smith et al. 2016 Smith, A. M.; Lee, A. A.; Perkin, S. The Electrostatic Screening Length in Concentrated Electrolytes Increases with Concentration. _J. Phys. Chem. Lett._ 2016, _7_ , 2157–2163
  * Fung and Perkin 2023 Fung, Y. K. C.; Perkin, S. Structure and anomalous underscreening in ethylammonium nitrate solutions confined between two mica surfaces. _Faraday Discuss._ 2023, _246_ , 370–386
  * Lee et al. 2017 Lee, A. A.; Perez-Martinez, C. S.; Smith, A. M.; Perkin, S. Underscreening in concentrated electrolytes. _Faraday Discuss._ 2017, _199_ , 239–259
  * Yuan et al. 2022 Yuan, H.; Deng, W.; Zhu, X.; Liu, G.; Craig, V. S. J. Colloidal Systems in Concentrated Electrolyte Solutions Exhibit Re-entrant Long-Range Electrostatic Interactions due to Underscreening. _Langmuir_ 2022, _38_ , 6164–6173
  * Attard 1993 Attard, P. Asymptotic analysis of primitive model electrolytes and the electrical double layer. _Phys. Rev. E_ 1993, _48_ , 3604–3621
  * Coupette et al. 2018 Coupette, F.; Lee, A. A.; Härtel, A. Screening Lengths in Ionic Fluids. _Phys. Rev. Lett._ 2018, _121_ , 075501
  * Härtel et al. 2023 Härtel, A.; Bültmann, M.; Coupette, F. Anomalous Underscreening in the Restricted Primitive Model. _Phys. Rev. Lett._ 2023, _130_ , 108202
  * Kumar et al. 2022 Kumar, S.; Cats, P.; Alotaibi, M. B.; Ayirala, S. C.; Yousef, A. A.; van Roij, R.; Siretanu, I.; Mugele, F. Absence of anomalous underscreening in highly concentrated aqueous electrolytes confined between smooth silica surfaces. _J. Colloid Interface Sci._ 2022, _622_ , 819–827
  * Rotenberg et al. 2018 Rotenberg, B.; Bernard, O.; Hansen, J.-P. Underscreening in ionic liquids: a first principles analysis. _J. Phys. Condens. Matter_ 2018, _30_ , 054005
  * Kjellander 2020 Kjellander, R. A multiple decay-length extension of the Debye-Hückel theory: to achieve high accuracy also for concentrated solutions and explain under-screening in dilute symmetric electrolytes. _Phys. Chem. Chem. Phys._ 2020, _22_ , 23952–23985
  * Coles et al. 2020 Coles, S. W.; Park, C.; Nikam, R.; Kanduc, M.; Dzubiella, J.; Rotenberg, B. Correlation Length in Concentrated Electrolytes: Insights from All-Atom Molecular Dynamics Simulations. _J. Phys. Chem. B_ 2020, _124_ , 1778–1786
  * Zeman et al. 2020 Zeman, J.; Kondrat, S.; Holm, C. Bulk ionic screening lengths from extremely large-scale molecular dynamics simulations. _Chem. Commun._ 2020, _56_ , 15635–15638
  * Cats et al. 2021 Cats, P.; Evans, R.; Härtel, A.; van Roij, R. Primitive model electrolytes in the near and far field: Decay lengths from DFT and simulations. _J. Chem. Phys._ 2021, _154_ , 124504
  * Ma et al. 2015 Ma, K.; Forsman, J.; Woodward, C. E. Influence of ion pairing in ionic liquids on electrical double layer structures and surface force using classical density functional approach. _The Journal of Chemical Physics_ 2015, _142_ , 174704
  * Hasted et al. 2004 Hasted, J. B.; Ritson, D. M.; Collie, C. H. Dielectric Properties of Aqueous Ionic Solutions. Parts I and II. _J. Chem. Phys._ 2004, _16_ , 1–21
  * de Souza et al. 2022 de Souza, J.; Kornyshev, A. A.; Bazant, M. Z. Polar liquids at charged interfaces: A dipolar shell theory. _J. Chem. Phys._ 2022, _156_
  * Conway and Marshall 1983 Conway, B.; Marshall, S. Some common problems concerning solvent polarization and dielectric behaviour at ions and electrode interfaces. _Aust. J. Chem._ 1983, _36_ , 2145–2161
  * Bonthuis et al. 2012 Bonthuis, D. J.; Gekle, S.; Netz, R. R. Profile of the Static Permittivity Tensor of Water at Interfaces: Consequences for Capacitance, Hydration Interaction and Ion Adsorption. _Langmuir_ 2012, _28_ , 7679–7694
  * Danielewicz-Ferchmin et al. 2013 Danielewicz-Ferchmin, I.; Banachowicz, E.; Ferchmin, A. Dielectric saturation in water as quantitative measure of formation of well-defined hydration shells of ions at various temperatures and pressures. Vapor–liquid equilibrium case. _J. Mol. Liq._ 2013, _187_ , 157–164
  * Adar et al. 2018 Adar, R. M.; Markovich, T.; Levy, A.; Orland, H.; Andelman, D. Dielectric constant of ionic solutions: Combined effects of correlations and excluded volume. _J. Chem. Phys._ 2018, _149_ , 054504
  * Underwood and Bourg 2022 Underwood, T. R.; Bourg, I. C. Dielectric Properties of Water in Charged Nanopores. _The Journal of Physical Chemistry B_ 2022, _126_ , 2688–2698, PMID: 35362980
  * Ben-Yaakov et al. 2011 Ben-Yaakov, D.; Andelman, D.; Podgornik, R. Dielectric decrement as a source of ion-specific effects. _J. Chem. Phys._ 2011, _134_ , 074705
  * Hubbard et al. 1979 Hubbard, J. B.; Colonomos, P.; Wolynes, P. G. Molecular theory of solvated ion dynamics. III. The kinetic dielectric decrement. _J. Chem. Phys._ 1979, _71_ , 2652–2661
  * Knight and Hub 2015 Knight, C. J.; Hub, J. S. WAXSiS: a web server for the calculation of SAXS/WAXS curves based on explicit-solvent molecular dynamics. _Nucleic Acids Research_ 2015, _43_ , W225–W230
  * Hansen et al. 2021 Hansen, J.; Uthayakumar, R.; Pedersen, J. S.; Egelhaaf, S. U.; Platten, F. Interactions in protein solutions close to liquid–liquid phase separation: ethanol reduces attractions via changes of the dielectric solution properties. _Physical Chemistry Chemical Physics_ 2021, _23_ , 22384–22394
  * Hansen et al. 2022 Hansen, J.; Pedersen, J. N.; Pedersen, J. S.; Egelhaaf, S. U.; Platten, F. Universal effective interactions of globular proteins close to liquid–liquid phase separation: Corresponding-states behavior reflected in the structure factor. _The Journal of Chemical Physics_ 2022, _156_
  * 38 Forsman, J.; Ribar, D.; Woodward, C. E. _submitted to Physical Chemistry Chemical Physics_
  * Torrie and Valleau 1980 Torrie, G. M.; Valleau, J. P. Electrical double layers. I. Monte Carlo study of a uniformly charged surface. _J. Chem. Phys._ 1980, _73_ , 5807–5816
  * Leote de Carvalho and Evans 1994 Leote de Carvalho, R.; Evans, R. The decay of correlations in ionic fluids. _Mol. Phys._ 1994, _83_ , 619–654