Consistent description of ion-specificity inbulk and at interfaces by solvent implicit
simulations and mean-field theory
Cite as: J. Chem. Phys. 153, 034103 (2020); https://doi.org/10.1063/5.0016103
Submitted: 02 June 2020 . Accepted: 22 June 2020 . Published Online: 15 July 2020
Alexandre P. dos Santos, Yuki Uematsu, Alexander Rathert, Philip Loche

J. Chem. Phys. 153, 034103 (2020); https://doi.org/10.1063/5.0016103© 2020 Author(s).

, and Roland R. Netz

153, 034103

The Journalof Chemical Physics

ARTICLE

scitation.org/journal/jcp

Consistent description of ion-specificity
in bulk and at interfaces by solvent implicit
simulations and mean-field theory
Cite as: J. Chem. Phys. 153, 034103 (2020); doi: 10.1063/5.0016103
Submitted: 2 June 2020 • Accepted: 22 June 2020 •Published Online: 15 July 2020
Alexandre P. dos Santos,1,2,a)and Roland R. Netz2,b)

Yuki Uematsu,3,b)Alexander Rathert,2,b) Philip Loche,2,b)

AFFILIATIONS
1Instituto de Física, Universidade Federal do Rio Grande do Sul, Caixa Postal 15051, 91501-970 Porto Alegre, RS, Brazil

2

Fachbereich Physik, Freie Universität Berlin, 14195 Berlin, Germany

3

Department of Physics, Kyushu University, 819-0395 Fukuoka, Japana)

Author to whom correspondence should be addressed: alexandre.pereira@ufrgs.br
Electronic addresses: uematsu.yuki@phys.kyushu-u.ac.jp; a.rathert@fu-berlin.de; ploche@physik.fu-berlin.de;
and rnetz@physik.fu-berlin.deb)

ABSTRACT

Solvent-implicit Monte Carlo (MC) simulations and mean-field theory are used to predict activity coefficients and excess interfacial tensions
for NaF, NaCl, NaI, KF, KCl, and KI solutions in good agreement with experimental data over the entire experimentally available concentration range. The effective ionic diameters of the solvent-implicit simulation model are obtained by fits to experimental activity coefficientdata. The experimental activity coefficients at high salt concentrations are only reproduced if the ion-specific concentration-dependent decrement of the dielectric constant is included. The dielectric-constant dependent contribution of the single-ion solvation free energy to theactivity coefficient is significant and is included. To account for the ion-specific excess interfacial tension of salt solutions, in addition to nonideal solution effects and the salt-concentration-dependent dielectric decrement, an ion-specific ion–interface interaction must be included.
This ion–interface interaction, which acts in addition to the dielectric image-charge repulsion, is modeled as a box potential, is considerably more long-ranged than the ion radius, and is repulsive for all ions considered except iodide, in agreement with previous findings andarguments. By comparing different models that include or exclude bulk non-ideal solution behavior, dielectric decrement effects, and ion–interface interaction potentials, we demonstrate how bulk and interfacial ion-specific effects couple and partially compensate each other. Our
MC simulations, which correctly include ionic correlations and interfacial dielectric image-charge repulsion, are used to determine effectiveion–surface interaction potentials that can be used in a modified Poisson–Boltzmann theory.
Published under license by AIP Publishing. https://doi.org/10.1063/5.0016103., s

I. INTRODUCTION
The ionic specificity observed in interfacial tension data,1–7 critical coagulation or micellization concentration data of colloidal andmicellar suspensions,8–16 surface interactions,17,18 and protein structure and kinetics19,20 still poses a number of questions.21–24 This isso because in all ion-specific phenomena that involve surfaces incontact with salt solutions, three distinct effects are intimately coupled: (i) Ions interact with interfaces not only by universal Coulomb
and image-charge forces but also by ion-specific force contributions that involve hydration effects and van der Waals interactions.

J. Chem. Phys. 153, 034103 (2020); doi: 10.1063/5.0016103
Published under license by AIP Publishing(ii) The bulk ion chemical potential depends on salt concentration
with pronounced variations among different ions, and this is at theheart of bulk ion specificity and sensitively influences the ionic surface adsorption equilibrium. (iii) Finally, the dielectric constant ofsalt solutions decreases with rising salt concentration in a manner that is again different for different ions. The dielectric constant
not only influences the bulk ion chemical potential but also theion–interface dielectric image-charge interactions. All three abovementioned effects have been amply treated separately in the literature, but an investigation of how these effects conspire to produce
ion-specific effects at interfaces has been largely missing.

153, 034103-1

The Journalof Chemical Physics

In early modeling approaches, the dispersion or van der Waals
interaction between ions and the air–water interface was suggested
to be responsible for ion-specific effects seen in air–water interfacial tension data.25 Ion-specific interactions between ions and theair–water interface that reproduced and explained the experimental trends were in water-explicit simulations first observed withpolarizable force fields,26,27 and such simulations also demonstrated
the enhanced interface activity of iodide.28–31 The contribution ofionic polarizability to ion–interface interactions was investigated ina number of subsequent theoretical and simulation works. As a general trend, it was seen that small ions with tightly bound hydrationwater are strongly repelled from the electrolyte–air interface, whilelarger ions with loosely bound hydration water are less strongly
repelled or can even adsorb at the interface.32–35
It was later shown that well-parametrized non-polarizable
force-field simulations can also account for ion-specific ion–interface interactions and, in conjunction with mean-field
Poisson–Boltzmann modeling, reproduce experimental excess surface tensions.36,37 More recently, the salt-concentration-dependent
air–water interfacial tension was directly extracted from waterexplicit molecular dynamics simulations and shown to reproduce
the experimental difference between acids, which lower the interfacial tension, and bases, which raise the interfacial tension.38 Whilewater-explicit simulations, in principle, account for hydration andnon-ideal solution effects, the salt-induced lowering of the dielectric
constant, and ion–surface interactions, provided that accurate force
fields are used, the in-depth understanding of ion-specific effects andthe description of experimental scenarios that go beyond idealized
planar interfaces require the use of coarse-grained simulations andanalytically tractable theories.
Besides simulations, there are several approximate methods
that can be used to predict electrolyte chemical potentials. Startingwith the pioneering work of Pitzer,39,40 various models were developed to describe the electrolyte activity by simple analytical expressions.41–43 The Mean Spherical Approximation (MSA) in fact very
accurately predicts the salt solution activity compared to simulations using the same underlying Hamiltonian.44–46 It is well known
that the dielectric constant of electrolyte solutions decreases withelectrolyte concentration,47 and consequently, it was shown that inorder to correctly predict the experimental activity of salt solutions
at high concentrations, this dielectric decrement must be taken into
account.48–53 Note that the dielectric decrement influences also thesingle-ion solvation free energy and thereby makes the solvation
free energy dependent on the salt concentration.52,54 There are several models that describe the dielectric constant of salt solutions:55
the heuristic Gavish and Promislow model56 describes experimental
data very accurately over the entire concentration range in terms of
two fit parameters, namely, the infinite-dilution excess ionic polarizability and the molten-salt dielectric constant at zero water content.
Using solvent-implicit simulations of interfacial ionic systems, it wasseen that non-ideal solution effects as well as dielectric interfacial
effects crucially influence the interfacial tension.57–59
In our work, we aim at integrating all these effects by the useof suitably constructed solvent-implicit simulation models and bymodified Poisson–Boltzmann (MPB) approaches for monovalent
electrolytes. As a simulation model, we use the so-called primitivemodel, where ions are described by hard charged spheres and water
is described by a medium with a homogeneous dielectric constant.

J. Chem. Phys. 153, 034103 (2020); doi: 10.1063/5.0016103
Published under license by AIP Publishing

ARTICLE

scitation.org/journal/jcp

TABLE I. Summary of the different MPB and MC models. For the MC models, thescreened image-charge OS potential is implicitly included.

MPB-1
MPB-2
MPB-3
MPB-4
MPB-5
MC-2
MC-3
MC-4
MC-5

OS

γ±

ε(cb )

No
Yes
Yes
Yes
Yes
...
...
...
...

Ideal
Ideal
Non-ideal
Ideal
Non-ideal
Ideal
Non-ideal
Ideal
Non-ideal

Constant
Constant
Constant
DependDependConstant
Constant
DependDependWe first compare experimental activity coefficient data for six different salts with grand-canonical Monte Carlo simulation (GCMC)
results and with MSA results over the entire concentration range.
The agreement between GCMC and MSA results is good for mostconcentrations, which is crucial because it allows us to use theclosed-form MSA chemical potential expression for our subsequentcalculation of interfacial tensions. The dielectric decrement of thesalt solutions is in the simulations, and MPB approaches accounted
for by the use of the Gavish and Promislow model, where the excess
ionic polarizabilities are taken from experimental data, while themolten-salt dielectric constant parameters are obtained by fitting theexperimental dielectric constant curves. The fit of the GCMC predictions to the experimental activity coefficient data yields in a second
step the effective ionic diameters.
Using the salt-concentration-dependent bulk dielectric constants and the effective ionic diameters determined by comparison
with experimental data, we then calculate excess surface tensions ofsalt solutions by means of Monte Carlo (MC) simulations that properly account for dielectric image-charge effects at the surface. The
comparison with experimental data demonstrates that an additional
ion-specific interaction between the air–water interface and the ionsis present, which is more long-ranged than the ion radius and whichis repulsive for all ions except iodide. Here, we also compare witha MPB model that accounts for image-charge effects [similar to theOnsager–Samaras (OS) theory60 ] but that also includes additional
ion-specific ion–interface interactions. We end with a discussionof the relative importance of non-ideal solution effects, the dielectric decrement of salt solutions, and ion–interface interactions that
include dielectric image-charge effects. Table I lists different MC and
MPB models that are compared to each other.
II. METHODS
A. Analytical model for activity coefficients
The mean ionic chemical potential, averaged over cations andanions, can be written as
μ± (cb ) = kB T ln (d3 cb γ± ),(1)

where cb is the bulk salt concentration in terms of the number density of ion pairs, kB T is the thermal energy, d is the mean ionic
diameter, and γ± is the mean ionic activity coefficient. Deviations ofγ± from unity indicate non-ideal solution effects due to interactions

153, 034103-2

The Journalof Chemical Physics

ARTICLEbetween ions and between ions and water.52 For a system of charged
hard spheres that interact via Coulomb interactions, the so-calledprimitive model, the mean ionic activity coefficient can be obtainedby adding two contributions, the ion–ion term, modeled with MSA
theory,44,48,57 and the ion–water contribution μiw ,√κd (1 + 2κd) − κd − (κd)2 8ϕ − 9ϕ2 + 3ϕ3ln γ± =+
+ βμiw , (2)24ϕ(1 − ϕ)3√
where κ = 8πq2 ℓB (cb )cb is the inverse Debye length;ℓB (cb ) =e2
4πε(cb )ε0 kB T

(3)

is the salt-concentration-dependent Bjerrum length; e is the elementary charge; q is the ion valency; ε(cb ) is the relative electrolyte
dielectric constant, which depends on the bulk salt concentrationand thereby makes the Bjerrum length and the Debye length saltconcentration dependent; ε0 is the dielectric permittivity of vacuum;and ϕ = cb πd3 /3 is the ionic volume fraction. The Bjerrum length atinfinite dilution cb = 0 is set to ℓB (0) = 0.71 nm, which is obtainedusing the pure water relative dielectric constant61 εw = 78.5 at room
temperature, 25 ○ C. The mean ionic activity coefficient depends on
the ion type via the ionic diameter d, which is different for different ion combinations, and via the salt-concentration-dependent
dielectric constant ε(cb ). The physical meaning of the mean ionic
diameter d is the minimal cation–anion separation at contact. Thisis so because the cation–anion interaction is attractive and therefore
dominant compared to cation–cation and anion–anion interactions.
The concentration-dependent ion–water contribution is obtained by
the Born solvation model,52
βμiw =

e2
1
1
(
− ),
4πε0 kB TaB ε(cb ) εw

(4)

where εw = ε(cb = 0) denotes the pure water dielectric constant bywhich the ion–water contribution vanishes in the infinite-dilution
limit. Here, aB is the ionic mean Born diameter, which is obtained
from experimental data for the solvation free energy in the infinitedilution limit according to62
ΔG =

1
e2
( − 1).
4πε0 aB εw(5)Note that the ionic mean Born diameter aB is not necessarily thesame as the mean ionic diameter d that appears in the ion–ion
interaction contribution to the activity coefficient. This reflects thefact that aB results from ion–water interactions and is modified
due to non-linear dielectric effects,63 while d reflects the closest
distance between cations and anions. Fitting results for all diameters are shown in Table II. In Appendix A, we compare our fitting
results with results for the NaCl solvation free energy obtained fromall-atom MD simulations.
B. Dielectric constant at finite salt concentration
The dielectric constant of salt solutions decreases significantlywith rising salt concentration in a fashion that depends on the salt
type. The experimental ion-specific dielectric decrement can be very

J. Chem. Phys. 153, 034103 (2020); doi: 10.1063/5.0016103
Published under license by AIP Publishingscitation.org/journal/jcp

well described by the heuristic Gavish–Promislow model,56
⎡
⎤
⎢
⎥
⎢
(εw − εms ) ⎥1⎢⎥⎥.
ε(cb ) = εw + (εw − εms )⎢−⎢3αp cb ⎥⎢ tanh ( 3αp cb )⎥⎢⎥εw − εms⎣⎦(6)

In the limit of low concentration cb , one obtains ε(cb ) ≃ εw + αp cb ,while in the limit of high concentration, one obtains ε(cb ) ≃ εms .
The model contains two fitting parameters: the dielectric constantof the molten salt, denoted as εms , and the excess ionic polarizability per ion pair αp in the infinitely dilute limit. Here, we take thevalues of αp from experimental measurements, listed in Table II,so that we are left with only a single fit parameter εms for each ion
pair.
C. Grand-canonical Monte Carlo simulations
for bulk salt solutions
To obtain the activity coefficient of a bulk salt solution as afunction of salt concentration, we perform grand-canonical simulations and thereby measure the concentration as a function of
the imposed chemical potential. The Born solvation part, given by
Eq. (4), is added to the chemical potential imposed in the simulations. Besides regular MC moves, where single ions are moved in
random directions, the GCMC method consists of insertions andremovals of ion pairs that conserve the system electroneutrality,
using the proper acceptance rules.64,65 In the present case, where
the dielectric constant depends on the salt concentration, we update
the dielectric constant, which enters the calculation of the systemenergy, every 10 MC steps according to the mean ion concentration. For high salt concentrations, where the dielectric constants are
low, we update the dielectric constant more frequently in order toaccelerate equilibration.
The Coulomb energy is calculated using Ewald summation
methods65 based on the energy expression
qj qk
⎧
N
N
⎪
⎪
⎪ℓB (cb ) ∑j=1 ∑k=j+1 r
jk
U/kB T = ⎨
⎪
⎪
⎪
⎩∞

for rjk ≥ d(7)

for rjk < d,

TABLE II. Summary of electrolyte parameters. The ionic diameters d are obtained byfits to the experimental activity coefficients in Fig. 3. The limiting dielectric constantsat zero water content, εms , are obtained by fits to the experimental dielectric data
in Fig. 2, while the excess ionic polarizabilities αp are extracted from experimental
data.47,78,81,82 The value of εms for NaF is taken to be the same as KF because of therestricted experimental concentration range of the dielectric data for NaF. The mean
ionic Born diameters aB are obtained from the Born model in Eq. (5) and experimental
data from Ref. 62.

Salt

d (nm)εms

αp (M− 1 )aB (nm)

NaCl
KCl
NaI
KI
NaF
KF

0.249
0.23
0.27
0.25
0.17
0.24

27.9
35
18.3
30.5
22.4
22.4

−11.59−10.02−14.58−14.69
−11.9
−12

0.369
0.406
0.392
0.433
0.340
0.371

153, 034103-3

The Journalof Chemical Physics

ARTICLE

scitation.org/journal/jcp

where rjk is the distance between ions j and k of valencies qj andqk , respectively, and N is the fluctuating total number of particles. Typical particle numbers are between 500 and 1000. The cubic
simulation box size is between Lb = 4 nm and Lb = 25 nm, dependingon the bulk ion concentration.
D. Canonical Monte Carlo simulations in slab
geometry
In the slab geometry, we perform canonical simulations at a
fixed total ion number and estimate the bulk ion concentration fromthe concentration at the mid plane between the two surfaces, separated by distance L. The simulation model is schematically shownin Fig. 1. Cations and anions have the same diameter d, which is
obtained from fits to experimental mean activity coefficient data
γ± , as will be explained in Sec. III. Ions interact with the interfaces
via a square-well box potential of strength αi if they approach theinterfaces closer than the distance h± = 0.5 nm. The interface interaction parameters αi are fitted by comparison with experimental
interfacial tension data. The interaction range h± = 0.5 nm, which
is significantly larger than the ionic radii that are in the range of d/2= 0.08 nm–0.13 nm, was previously shown to be in qualitative agreement with interface-ion interaction profiles extracted from waterexplicit molecular dynamics simulations and furthermore allows us
to robustly fit experimental interfacial tension data.66 In one test
simulation, we show interfacial tension results for the choice h±= d/2 [the green line in Fig. 4(a)] and demonstrate that this interaction range is too small and cannot be used to describe the experimental data.
The Coulomb energy for MC slab simulations includes direct
ion–ion as well as ion–interface image-charge interactions. There
are various methods for taking into account surface polarization
effects due to the difference of the dielectric constants of the salt
solution and the surrounding air.67–72 We choose the method developed in Ref. 73, which extends the modified Ewald sum and considers an infinite number of image charges. We therefore use the energy
expression
qj qk
⎧
⎪
ℓB (cb ) ∑Nj=1 ∑Nk=j+1
⎪
⎪
⎪
rjk
⎪
⎪
⎪
l
⎪
⎪
χ(c
)
ℓ
(c
)
B
b
b
⎪
∞
N
N qj qk⎪⎪+ ∑l=1∑j=1 ∑k=1 ′
⎪
⎪
2
rjk
⎪
⎪
⎪
⎪
l
U/kB T = ⎨
q
χ(c
)
ℓ
(c
)
j qk
B b
b
N
N
⎪
+ ∑∞
∑j=1 ∑k=1 ′′
⎪
l=1
⎪
⎪
2
r
⎪
jk
⎪
⎪
⎪
⎪
⎪
+
U
/k
T
for
r
≥
d
s
B
⎪
jk
⎪
⎪
⎪
⎪
⎪
⎪
⎪
⎩∞ for rjk < d, zj , < 0 or zj > L,

(8)

where χ(cb ) = [ε(cb ) − 1]/[ε(cb ) + 1] is the dielectric contrast between
the slab interior, which is filled with electrolyte, and the surrounding
air. The surface energy U s is given byN

Us /kB T = ∑ αi (j)[1 − θ(zj − h± ) + θ(zj − L + h± )].

(9)

j=1

The function θ(z) is the Heaviside function, and the surface interaction strength αi (j) for the jth ion depends on ion type i. N = 1000is the total number of particles (500 cations and 500 anions), rjk is

J. Chem. Phys. 153, 034103 (2020); doi: 10.1063/5.0016103
Published under license by AIP Publishing

FIG. 1. Schematic geometry of the model for ions near air–electrolyte interfaces.
(a) In the electrolyte medium, the dielectric constant depends on the bulk ionic concentration, ε(cb ), while in air, the dielectric constant is set to 1. (b) The interface box
potential accounts for ion–interface interactions that are due to interfacial hydration effects. (c) Exemplary ionic density profile for NaI at a bulk salt concentrationof cb = 0.53 M from MC simulations using the MC-5 model.

the distance between ions j and k, and rjk′ and rjk′′ are the distances
between charge j and image charges k located in regions z < 0 andz > L, respectively. The simulation box widths are Lx = Ly in
the lateral directions (depending on the ionic concentration) and
Lz = 20Lx in the direction normal to the interfaces (see Refs. 73 and74 for more details). The interface separation is chosen as L = 20κ−1 ,where κ is the inverse Debye length, in order to prevent interactions between the two interfaces. Thus, the Coulomb energy used in
the MC simulations takes into account the dependency of the Bjerrum length and the electrolyte dielectric constant on the bulk ion
concentration.

153, 034103-4

The Journalof Chemical Physics

ARTICLEThe excess interfacial tension is determined by the Gibbs
adsorption isotherm
Δσ = −2 ∫ Γ± dμ± ,(10)where Γ± is the mean ionic excess per unit area, which follows fromthe ionic concentration profiles obtained in the MC simulations asN− cb L/2,
4L2x

Γ± =

(11)

E. Modified Poisson–Boltzmann theory
The MPB theory provides an analytical route to predict theinterfacial tension of electrolytes including the effects of interface–ion interactions and dielectric decrement effects. For an interface
located at z = 0 and the electrolyte present at z > 0, the short-ranged
effective interaction potential U i (z) between an ion of type i and theinterface can be written as(12)

where αi is the ion-specific amplitude of the box potential, as in
the MC simulations, and the first term reflects screened imagecharge interactions.75,76 The mean ion diameter d and the box potential range h± are the same as in the MC simulations. The electrolyte dielectric constant depends on the ionic bulk concentration
cb according to Eq. (6), and thereby, χ, κ, and ℓB depend on cb as well.
In the limit when the dielectric constant is salt independent and theionic diameter d is taken to be zero, the first term in Eq. (12) reduces
to the interface potential used in the original Onsager–Samaras (OS)theory.60
We include the ion–interface interaction potential into the
MPB equation asε(cb )ε0d2 ψ= −ecb ∑ qi e−eqi ψ(z)/kB T−Ui (z)/kB T .
dz2
i

(13)

Here, we neglect any dependence of ε on the local ion concentrationand consider it to be uniform.
The individual ionic surface excesses follow as∞

Γi = cb ∫0(e−eqi ψ(z)/kB T−Ui (z)/kB T − 1)dz.

(14)

Charge neutrality dictates Γ+ = Γ− ; we thus define the mean ionicsurface excess Γ± = Γ+ = Γ− . The Gibbs adsorption isotherm[Eq. (10)] can be written in terms of the activity coefficient as
Δσ = −2kB T ∫cb0Γ± (1 + cb

d ln γ± dcb
)
.
dcb
cb

J. Chem. Phys. 153, 034103 (2020); doi: 10.1063/5.0016103
Published under license by AIP Publishing

In the ideal limit, one has cb d ln γ± /dcb = 0, while in the general nonideal case, we obtain from Eq. (2),√d ln γ± κd 1 + 2κd − κd(1 + κd) 2ϕ(4 − ϕ)
=
cb
+
√
dcb
(1 − ϕ)4
16πcb d3 1 + 2κd
1 √
κd
κd dε
−
− 1 − 2κd]
[ 1 + 2κd + √8πd3
2ε dcb
1 + 2κd
2
cb dε
e
,
(16)
−
4πε0 kB TaB ε2 dcb
where we obtain from Eq. (6) the result

where cb is the mean ionic bulk concentration far from the interfaces
in the middle of the simulation box. Using for the non-ideal chemical potential Eq. (1), Eq. (10) can be straightforwardly integrated
numerically.

Ui (z)
eκd ℓB (cb ) −2κz= χ(cb )e+ αi θ(h± − z),kB T
1 + κd 4z

scitation.org/journal/jcp(15)(εw − εms )2
dε
=
−
dcb
3αp cb23αp.
2
3αp cb[sinh()]
εw − εms(17)

To discuss the relative importance of non-ideal solution, ioncorrelation, and salt-concentration-dependent dielectric effects, weconsider five different MPB models and four different MC models, in which we separately include or neglect the non-ideal chemical potential, the concentration dependence of ε, and, for the MPB
model only, the OS potential [the first term in Eq. (12)]. The different models are summarized in Table I. For the MPB-1 model, an
analytic solution is available.66 The MPB-2 model is identical to theOnsager–Samaras theory in the limit of κd → 0 and αi → 0, in which
limit the surface tension can be exactly calculated as60
Δσ = −

2πkB T ∞ nn hn+1
(ln h + gn ),∑ℓB 2 n=1 n!n! n + 1(18)

where gn = ln n + 1 + 2γ − 1/(n + 1) − ∑nm=1 (1/m), h = χκℓB /2, andγ is the Euler constant. A comparison of the Onsager–Samaras theory with MPB modeling results is shown in Fig. 10 in Appendix B.
For models MPB-2 to MPB-5, we solve Eq. (13) by the shooting
method and then integrate Eq. (15) using the trapezoidal rule withintegration steps of Δcb = 10 mM.
III. RESULTS
A. Bulk effects
In Fig. 2, we compare experimental dielectric constant datawith the heuristic model [Eq. (6)]. The values for the excess ionicpolarizabilities, αp , are directly taken from experiments, where weaverage the values given in Refs. 47, 78, 81, and 82 (see Table II).
The resultant fit values for εms are given in Table II. The fit of εms
for NaF is not conclusive because the maximum solubility is around
1 M, and thus, the available data range is very small; we therefore use
for NaF the same εms value as determined for KF.
In Fig. 3, we compare experimentally measured activity coefficients for various salts83 (filled symbols) with GCMC simulation
results (open symbols) as well as with MSA theory (solid and broken lines). In both theoretical models, we include the Born solvation
contribution βμiw from Eq. (4). In Fig. 3(a), we show results forthe sodium salts, and in Fig. 3(b), we present results for the potassium salts. Note the different ordering of the salts with respect to theanions: in Fig. 3(a), the expected ordering iodide–chloride–fluoride(from top to bottom) is obtained, while in Fig. 3(b), the ordering

153, 034103-5

The Journalof Chemical Physics

ARTICLE

scitation.org/journal/jcp

FIG. 2. Dielectric constants as a function of salt concentration. Symbols are experimental data,55,77–80 and lines correspond to the Gavish–Promislow model [Eq. (6)] withparameters determined by a non-linear fitting procedure (see Table II). (a) Data for NaF, NaCl, and NaI. (b) Data for KF, KCl, and KI.

is partially reversed. Such ion-specific series reversals are well documented in the literature,52,86 have their origin in nonlinear ionic
hydration effects, and occur for all different kinds of experimental
observables.87
The only remaining fitting parameter in the modeling is theionic diameter d, which we determine by a fit of the GCMC results tothe experimental data; for the MSA predictions, we use the same values of d. It is seen that the MSA model agrees well with GCMC andreproduces the experimental activity coefficients in a satisfactory
fashion over the entire available concentration range. The brokenlines show MSA theory predictions for the case where the dielectric
constant is fixed at the pure water value and the ionic diameters are
the same as used for the solid lines. The difference between these
two MSA predictions is immense, which illustrates the importanceof taking the salt-induced modifications of the dielectric constantinto account.

FIG. 3. Mean ionic activity coefficient as a function of salt concentration for (a) NaF, NaCl, and NaI, and (b) KF, KCl, and KI. The experimental data are represented by fullsymbols.83 The empty symbols represent the GCMC simulations combined with the Born solvation model including the salt-concentration-dependent dielectric function, fromwhich the ion diameters d are obtained by fits to the experimental data. The solid lines represent MSA theory [Eq. (2)] using the same ionic diameters, dielectric molten-saltparameters εms , and Born solvation contribution. The dashed lines represent MSA theory for the same ion diameters d but for fixed dielectric constant corresponding to purewater. The units were transformed from molality to molarity using experimental solution mass density data for KCl84 and NaCl.85 As the correction is negligible for KCl andNaCl, we assume the same to be true also for the other salts (for which no data are available).

J. Chem. Phys. 153, 034103 (2020); doi: 10.1063/5.0016103
Published under license by AIP Publishing

153, 034103-6

The Journalof Chemical Physics

The parameters obtained for all salts are summarized inTable II. The agreement between simulation results and MSA theory
is satisfactory; we will therefore use the MSA result for the activity coefficient [Eq. (2)] when needed later to calculate the tensionof the electrolyte–air interface according to Eq. (10). The deviations
observed between MSA and simulations for higher concentrations,as seen for NaI in Fig. 3(a), can be explained by the low dielectric
constant of NaI solutions at high concentrations, which leads to ahigh electrostatic coupling constant, as is explained in Appendix C.
Interestingly, we note that when one neglects the ion solvation
contribution to the activity coefficient, the agreement between
experimental data and model prediction looks better, while the fitted diameters increase significantly, as discussed in Appendix D(Fig. 12).
The fit results for the ionic diameter d demonstrate that there
is no unique way to split the ionic diameters into individual ionicradii. Consider, for example, the fitted diameters d for NaCl andKCl (see Table II). For the salts involving Cl– , the Na+ radius would
be higher than the K+ radius, while for the salts involving F– , theopposite trend is obtained and the K+ radius is larger than the Na+radius. As mentioned before, this ion specific reversal of effective
ionic radii is well-known in the literature.52,88 The large difference
between the ionic Born diameters aB and the interaction ionic diameters d as extracted from fits of the experimental activity coefficients

ARTICLE

scitation.org/journal/jcp

is in Appendix A discussed in the context of non-linear dielectriceffects, which are extracted from water-explicit all-atom molecular
dynamics simulations. Since the solvation even of monovalent ionsinvolves considerable non-linear dielectric contributions, the Born
diameter aB has only limited physical meaning.
B. Interfacial effects
We now turn to interfacial effects. We calculate the excess surface tension by a combination of MC simulations and MSA theory,
where we integrate Eq. (10) based on ion density profiles from MC
simulations, using for the non-ideal chemical potential the MSA
expression [Eq. (1)]. We first consider the MC-5 model that includes
non-ideal bulk as well as concentration-dependent dielectric effects
using parameters that were obtained from fits to experimental activity coefficients (see Table II). We first neglect ion–interface interactions and set αNa = αCl = αI = 0. The MC simulation results
are represented by red solid lines in Figs. 4(a) and 4(b) and compared with experimental data from different sources. The simulation
results significantly underestimate the experimental data for NaClin (a). The green line in Fig. 4(a) shows simulation results wherea steric (i.e., infinitely high) ion–surface repulsion is considered,and we use d/2 = 0.1245 nm as the distance of closest approachfor Na+ as well as Cl− ions. It is seen that steric repulsion from the

FIG. 4. Comparison of experimental and MC simulation excess interfacial tension data. (a) NaCl: squares, right triangles, and diamonds represent the experimental dataof Matubayasi et al.,6 Jones and Ray,2 and Jarvis and Scheiman,3 respectively. (b) NaI: squares, diamonds, and circles represent the experimental data of Matubayasi
et al.,6 Jarvis and Scheiman,3 and Weissenborn and Pugh,4 respectively. (c) KCl: squares, right triangles, and circles represent the experimental data of Matubayasi et al.,6
Jones and Ray,2 and Weissenborn and Pugh,4 respectively. (d) KI: down triangles and diamonds represent the experimental data of Ali et al.89 and Jarvis and Scheiman,3
respectively. Solid lines denote MC simulation results using the MC-5 model. The dashed line in (c) represents MC simulation results using the MC-3 model, where thedielectric constant equals pure water while a non-ideal chemical potential is used. The dotted-dashed line in (c) denotes simulations using the MC-4 model where an ideal
chemical potential is used, but a concentration-dependent dielectric constant is employed.

J. Chem. Phys. 153, 034103 (2020); doi: 10.1063/5.0016103
Published under license by AIP Publishing

153, 034103-7

The Journalof Chemical Physics

ARTICLE

scitation.org/journal/jcp

TABLE III. Summary of the interfacial affinity parameters αi obtained from the fits in
Figs. 4 and 6.

MC-5
MPB-1
MPB-5

FIG. 5. Comparison of experimental interfacial tension data for NaF4 and MC
simulation results using the MC-5 model for all salt types considered in this work.

interface with a range given by the ionic radius d/2 is not sufficient
to correctly account for the NaCl experimental data. The violet curve
represents simulation results where we use an ion–interface interaction range of h± = 0.5 nm (as explained in Sec. II) and obtain byfitting αCl = 0.7 while fixing αNa = 0. Alternatively, the blue curveis the best fit when we enforce αNa = αCl , in which case we obtain
αNa = 0.35 and αCl = 0.35. It is seen that both parameter sets
reproduce very well the experimental data. In consistency with ourprevious modeling results66 and with our previous molecular
dynamics simulation results,37 we choose the latter parameter set
where sodium and chloride are equally repelled from the interface
for our further calculations. Note that this repulsion acts on top of
the dielectric image-charge repulsion, which is included explicitly in
the MC simulations.
By fixing αNa = 0.35, we obtain αI = −0.4 by a fit of our simulation results (purple curve) to the experimental data of NaI in
Fig. 4(b). Here, it is seen that choosing αNa = 0.35 and αI = 0.35
(blue curve) overestimates the experimental data, while choosing

Na

K

Cl

F

I

0.35
1.2
0.8

0.35
0.9
0.8

0.35
1.2
0.8

0.35
0.9
0.9

−0.4
0.3
−0.4

αNa = 0 and αI = 0 (red curve) describes the experimental data quite
well. We thus conclude that the iodide ion is slightly attracted to theair–electrolyte interface with a strength that amounts to roughly athird of the thermal energy, in agreement with previous simulation
reports28,31 and experiments.90 For the KCl and KI data in Figs. 4(c)and 4(d), we fix αCl = 0.35 and αI = −0.4 and choose αK = 0.35, which
gives a reasonably good comparison with the experimental data. We
note that because of the considerable scattering between the experimental data from different sources, a fit of the experimental data is
not meaningful.
In Fig. 5, we compare experimental data for NaF with MC
results for which we choose αF = 0.35, and here, a fit is not meaningful as the experimental data extend only to 1 M. In Fig. 5, we
also include a summary of the MC results for all different salts
considered in this paper, including KF (for which no experimental data are available) using the previously obtained parametersαK = 0.35 and αF = 0.35. The difference in the excess interfacial tension data between different salts is significant at high concentrations,and we note that these differences are caused by the combinationof ion–interfacial interactions as well as dielectric and non-ideal
effects. We conclude this section by noting that sodium, potassium, fluoride, and chloride are equally repelled from the interface, while the larger iodide ion is slightly attracted to the interface.
These ion–interface interactions do not primarily reflect dielectric image effects, since these are included in the MC simulations,rather they stem from ion-specific interfacial hydration effects that

FIG. 6. Comparison of the excess surface tension between MC simulations using the MC-5 model and two different MPB model results. Solid lines correspond to MC
simulation results, whereas dashed lines correspond to MPB results. (a) MPB-1 model results are shown that neglect the image-charge OS potential, the non-ideal activity,
and the concentration-dependent dielectric constant. (b) MPB-5 model results are shown that include the OS potential, the non-ideal activity, and the concentration-depending
dielectric constant. The fitted affinity parameters αi of the two different MPB models are summarized in Table III. For all models, we use εw = 78.5, T = 298 K, andh± = 0.5 nm as well as d, αp , and εms , as obtained in Sec. II A.

J. Chem. Phys. 153, 034103 (2020); doi: 10.1063/5.0016103
Published under license by AIP Publishing

153, 034103-8

The Journalof Chemical Physics

go beyond the solvent-implicit point-like dielectric ion model. The
important message here is that these effects can be accurately
accounted for by adding short-ranged ion–interface interaction
potentials into the MC simulations with interaction strength parameters that are obtained from fits to experimental interfacial tension
data.
C. Comparison of MC simulations and modified
PB models
So far, we compared MC simulations combined with MSA theory to experimental data in bulk and at interfaces. A very successful
and standard model to predict the behavior of ions at interfaces is
the Poisson–Boltzmann theory, which neglects correlations between
ions but otherwise treats the long-ranged electrostatic interactions
between ions exactly and has been very successfully applied to allkinds of different systems.91 At planar interfaces, the main shortcoming of PB approaches is that they neglect ion fluctuation effects
but also the image-charge repulsion of ions from low-dielectric surfaces. This is so because the PB approach treats a smeared-out homogeneous ion cloud that does not produce an electric field that penetrates into the substrate. We here modify the PB approach such as
to include the dielectric and other interfacial effects explicitly andcompare this MPB model with our MC simulations. As explained inSec. II, the MPB model includes image-charge repulsion via a heuristically added screened image-charge potential and is similar in spirit
to the original Onsager–Samaras theory.
In Fig. 6, we compare MPB theory (broken lines) with MC
simulation results (solid lines): in Fig. 6(a), we use the version
MPB-1 that neglects the dielectric OS potential, non-ideal solution,and dielectric decrement effects, and in Fig. 6(b), we use the version MPB-5 that accounts for all these effects. Note that the MCdata come from model MC-5 that includes non-ideal solution anddielectric decrement effects.
Our fitting strategy of the ion–interface interaction strengths αiin the two MPB models is the same as used previously when we fitted the MC simulations to the experimental data. First, we fit αNa

ARTICLE

scitation.org/journal/jcp

and αCl for αNa = αCl to the MC data for NaCl (red line in Fig. 6).
Then, we fit αK by comparison with the MC data for KCl (dark green
line in Fig. 6), and after that, we fit αI by comparison with the MCdata for NaI (purple line in Fig. 6). Finally, we fit αF by comparison
with the MC data for NaF (blue line in Fig. 6). All other parameters such as d, αp , and εms are taken as obtained in Sec. II A. We
find fair agreement between the MC and MPB-1 results in Fig. 6(a),and the agreement with the MPB-5 results in Fig. 6(b) is in fact very
good.
The resulting fit values for αi are summarized in Table III. We
observe that the MPB-1 αi values are significantly larger than theMC-5 values, which is easily explained by the fact that the MPB1 model does not include image-charge repulsion effects. To makeup for this, the image-charge repulsion is accounted for by the box
potential, which increases the interaction parameter. Note that theagreement between the MPB-1 model and the MC-5 simulation predictions for the interfacial tension in Fig. 6(a) is not bad, whichindicates that the MPB-1 model, which is analytically solvable, isquite useful and accurate in practice. The αi fit values for the MPB5 model are significantly reduced compared to the MPB-1 model,
which is understandable because now the image-charge repulsion is,according to Eq. (12), taken into account explicitly, but we still notea significant difference from the MC-5 fit results. We attribute thereason for this to the fact that the screened image-charge potential in
Eq. (12) underestimates the actual strength of the ionic image-charge
repulsion from the interface due to the effects of ion correlations. For
NaCl solutions, we obtain αNa = αCl = 1.2 using the MPB-1 model,
which is slightly larger than the fit values obtained previously using
the same model by comparison with molecular dynamics simulation results, in which case αNa = 1.16 and αCl = 0.98 were obtained.66
We note that even though we use the same αi values for NaCl and
KCl in the MPB-5 model, the predicted surface tension curves differ
from each other, similar to the MC-5 results. The reason for this isthat the values for the ion diameter d, the Born diameter aB , and thedielectric parameters αp and εms for NaCl and KCl slightly differ, assummarized in Table II.

FIG. 7. Decomposition of the effectsof the image-charge OS potential, thenon-ideal activity, and the concentrationdependent dielectric constant on theinterfacial tension. (a) Comparison of all
five MPB models. We use αK = αCl= 0.8 with corresponding d, αp , and εms
as appropriate for KCl. (b) Comparisonof all four MC models. We use αK = αCl= 0.35 with corresponding d, αp , and εms
as appropriate for KCl.

J. Chem. Phys. 153, 034103 (2020); doi: 10.1063/5.0016103
Published under license by AIP Publishing

153, 034103-9

The Journalof Chemical Physics

FIG. 8. Ion concentration profiles for a 0.53 M KI solution. The blue and red solidlines are the K+ and I− concentrations calculated using the MC-5 model, whereas
the blue and red dashed lines are calculated using the MPB-5 model.

In the following, we decompose the effects of the OS screened
image-charge potential, the non-ideal bulk activity, and the saltconcentration dependence of the dielectric constant on the interfacial tension, and we do this separately for the MC simulations andfor the MPB model. Figure 7(a) shows the interfacial tension calculated from all five different MPB models for fixed αK = αCl = 0.8. For
the MPB-1 model, the result (orange solid line) is strictly linear in thesalt concentration and is below all other MPB model results because
MPB-1 does not include the ion–interface dielectric repulsion asmodeled by the OS potential. The results of the MPB-2 and MPB4 models, which both include the OS potential but neglect non-ideal
solution effects, are almost the same, meaning that the bare effect of
the concentration-dependent dielectric constant (i.e., in the absenceof non-ideal solution effects) is minor. The result of MPB-3 has alower slope than the result of MPB-2, demonstrating the effect of
the non-ideal chemical potential. However, this decrease in slope
is partly canceled when including in addition the concentrationdependent dielectric constant, as demonstrated by comparison of
the MPB-5 result with the MPB-3 result. We conclude that ion–interface dielectric repulsion increases the interfacial tension, while
non-ideal solution effects and salt-concentration-dependent dielectric effects partly cancel each other. The comparison of the four
different MC models with αK = αCl = 0.35 in Fig. 7(b) exhibits very
similar trends.
In Fig. 8, we compare the ion density profiles obtained from
MC-5 simulations and the MPB-5 model. We show profiles for a KI
solution of bulk salt concentration cb = 0.53 M, where in the simulations the bulk salt concentration is determined from the concentration values far away from the interfaces. The MPB-5 model profiles
agree quite accurately with the MC simulation profiles, with somedifferences seen for the K+ profile near z = h± , where the box potential produces a discontinuity, which is more drastic for the MPB-5
model since the amplitude αK is higher, as shown in Table III. We
conclude that the MPB-5 model not only predicts interfacial tension data accurately but also gives a good account of the ion density
profiles.
IV. CONCLUSIONSThe excess air–water interfacial tensions of the six salt solutions NaF, NaCl, NaI, KF, KCl, and KI are calculated by means

J. Chem. Phys. 153, 034103 (2020); doi: 10.1063/5.0016103
Published under license by AIP Publishing

ARTICLE

scitation.org/journal/jcp

of solvent-implicit MC simulations that model the ions as charged
hard spheres. The MC simulations include correlation and fluctuation effects and account for dielectric image-charge effects thatstem from the dielectric contrast between air and aqueous salt solutions. The effective ionic interaction diameters d are obtained by fitsof simulated salt solution activity coefficients to experimental data.
This comparison also demonstrates that the dielectric decrementof salt solutions crucially influences the activities at high concentration and can be quantitatively described by a heuristic modelwith only two parameters that are extracted from experimental
data.
In solvent-implicit models, the ion radius is an effective parameter that accounts for ion-specific bulk hydration effects. This effective ion radius is not directly related to the ion crystal radius. Thisis most vividly demonstrated by the ion specific reversal seen in theexperimental bulk activity coefficient data in Fig. 3: in conjunction
with sodium as a counterion, the ion radius increases from fluoride
over chloride to iodide, as one would expect; in conjunction withpotassium as a counterion, however, the ion radius decreases fromfluoride to chloride.
By comparison of simulated and experimental interfacial tension data, we demonstrate that besides the dielectric decrement,
non-ideal solution effects, and image-charge interfacial repulsion,
we need to include an ion-specific ion–interface interaction, which
has a range that is significantly larger than the ion diameter andwhich we fix at 0.5 nm, guided by ion–interface interaction profiles
previously extracted from water-explicit molecular dynamics simulations. The ion–interface interaction parameters for Na+ , K+ , Cl– ,and F– are repulsive and given by α = 0.35 in units of kB T, whilefor I– , the interaction is attractive and given by α = −0.4. In theinterfacial MC simulations, the bulk chemical potential is needed for
calculating the interfacial tension from the interfacial ionic excesses.
For this, the chemical potential expressions from MSA theory areused, which are validated by bulk MC simulations.
In order to construct an even more coarse-grained description of ion-specific interfacial effects, we construct a mean-field
MPB model that besides non-ideal solution and dielectric decrement
effects also includes ion–interface interactions. The dielectric imagecharge repulsion of an ion from the low-dielectric air half space,which is missed in mean-field models at planar interfaces, is heuristically added to the ion potential by a screened image-charge potential,
very much in the spirit of the original Onsager–Samaras theory. The
strength of the ion–interface interaction potential is then obtainedby comparison of the MPB model with the MC simulation results.
MC and MPB predictions for interfacial tensions and ionic density
profiles agree quite nicely, while the effective ion–interface interaction in the MPB models tends to be more repulsive by not morethan 0.5kB T, which presumably comes from the neglect of correlation effects in MPB theory. By separately including and excluding
the effects of non-ideal electrolyte activity and the concentrationdependent dielectric decrement, we show that these two effects are
not additive and partially cancel each other.
Our MPB theory is a versatile, cheap, and quite accurate way of modeling ion-specific effects at interfaces, which
accounts for dielectric image-charge effects by a heuristic potential. The comparison with MC simulation results demonstrates
that the neglect of correlation effects, which is inherent to meanfield theories, is compensated by a slight shift of the effective

153, 034103-10

The Journalof Chemical Physics

ARTICLE

scitation.org/journal/jcp

ion–interface interaction potential strength. This explains why such
MPB models with suitably adjusted effective ion–surface interaction
potentials are able to describe experimental features very well.37,87
The minimal MPB model to consistently describe ion-specific effects
in bulk and at interfaces thus needs to include the following features: ion-type dependent chemical potentials (which are accurately
described by the closed-form MSA expression), ion-type dependent dielectric decrements (which are accurately described by theGavish–Promislow model), as well as short ranged and ion-type
dependent ion–interface potentials, with strengths that are adjusted
so as to account for dielectric image-charge effects as well as forinterfacial hydration effects.
ACKNOWLEDGMENTS
This work was supported by the Brazilian Research Council(CNPq) via Grant No. 302720/2018-9, FAPERGS, the Alexandervon Humboldt Foundation through the CAPES program, the MaxPlanck Water Initiative, and the Deutsche Forschungsgemeinschaft(DFG) via Grant No. NE 810/11. Y.U. was supported by JSPS KAKENHI, Grant No. 20K14430. R.R.N. acknowledges support from theInfosys Foundation.
APPENDIX A: ALL ATOM MOLECULAR DYNAMICSSIMULATIONS WITH EXPLICIT WATER: SINGLE-ION
SOLVATION FREE ENERGIES
Here, we calculate the solvation free energy F of single Na+
and Cl– ions by thermodynamic integration using all-atom molecular dynamics simulations. The ion Lennard-Jones (LJ) parameters
are taken from Ref. 92, and the used water model is SPC/E. The LJ
and Coulomb contributions to the free energy are calculated separately so that F = F LJ + F Coul . The explicit procedure and simulation
parameters are explained in Ref. 63. The Coulomb part of the freeenergy does not only contain the linear polarization contribution
as described by the Born solvation model given in Eq. (5), which
is proportional to the square of the ion charge, but also a term linear in the ion charge, which corresponds to an electrostatic potential
contribution, and non-linear dielectric terms. All these terms can be
described by the non-linear dielectric model according toFCoul = Φq + Aq2 + Bq3 + Cq4lin= Φq + FCoul+ Bq3 + Cq4 .

(A1)

Herein, A denotes the coefficient of the linear dielectric contribution, Φ is the electrostatic potential inside the ionic cavity due towater orientation effects, and B and C are non-linear coefficients.
To obtain the four coefficients, we calculate the free energy for different fractions of the full ion charge and fit the resulting Coulomb
free energy to Eq. (A1). We first fit up to second order using data
within a range of ±0.2e and afterward fit the remaining coefficientswith all the data. The fits are shown in Fig. 9, and Table IV shows
the fitted coefficients, which for Cl– agree with our previous work.63
The errors are estimated by varying the range from ±0.1e to ±0.3e
and taking the standard deviation. From the results of the fit, we canestimate the ionic Born diameters aB by using the Born equation, as
given in Eq. (5), together with the SPC/E dielectric constant εSPC/E= 71. Table V shows the solvation free energies F and theJ. Chem. Phys. 153, 034103 (2020); doi: 10.1063/5.0016103
Published under license by AIP Publishing

FIG. 9. Coulomb contribution to the free energy for different fractions of a fullcharged ion. The continuous line shows the fit according to Eq. (A1).

corresponding resultant values of aB separately for the Na+ and Cl–
ions. The two diameters are quite similar to each other and agree
reasonably well with the fit to experimental data in Table II. However, the Born model only describes the linear part of the Coulomb
part of the solvation free energy. When we fit only the Coulombic
part F Coul by the Born model, we obtain the modified ion diameters
aCoul , which differ only slightly from the ionic Born diameters aB .
lin
However, when we fit only the linear Coulombic part FCoul= Aq2 by
TABLE IV. Coefficients according to Eq. (A1) in units of kB T.

Ion

Φe

Ae2

Be3

Ce4

Na+
Cl–

19 ± 214.8 ± 0.4−181 ± 6−109.0 ± 0.8

62 ± 318.6 ± 0.5−48 ± 7
−11.5 ± 1

TABLE V. Solvation free energies (kB T) and corresponding ionic diameters (nm)according to Eq. (5).

Ion

F

aB

F Coul

aCoul

lin
FCoul

alin
Coul

Na+
Cl–

−146.0−144.2

0.376
0.381

−148.8−152.8

0.369
0.360

−181
−109.0

0.304
0.504

153, 034103-11

The Journalof Chemical Physics

ARTICLE

scitation.org/journal/jcp

as extracted from a fit of the experimental activity coefficients istherefore not surprising.
APPENDIX B: COMPARISON OF MODIFIED PB MODELWITH OS THEORY
Figure 10 shows a comparison of the MPB-2 model results with
the Onsager–Samaras theory. The red line is the Onsager–Samaras
theory [Eq. (18)], which is equal to the MPB-2 model with d = 0 nm
and αK = αCl = 0 (the dotted line). With d = 0.230 nm (the brokenline), the MPB-2 model shows a small difference from the Onsager–
Samaras theory, but still the experimental data significantly deviate
from the broken line. The black solid line is the fit of the MPB-2
model with the experimental data where d = 0.230 nm and αK = αCl= 0.8. Here, good agreement with the experimental data is observed.
APPENDIX C: VALIDITY OF MSA THEORY

FIG. 10. Comparison of the MPB-2 model results with the Onsager–Samaras theory. The black solid, broken, and dotted lines follow from the MPB-2 model fordifferent values of d, αK , and αCl as described in the legend. The red solid line
corresponds to the Onsager–Samaras theory calculated by evaluating the infinite
series in Eq. (18) up to n = 100, which perfectly agrees with the MPB-2 model
results for d = 0 nm and αK = αCl = 0 (black dotted line). The experimental data
deviate significantly from the Onsager–Samaras theory and can be described by
the MPB-2 model results when using repulsive ion–interface interactions with astrength of αK = αCl = 0.8 (solid line).

the Born model, we obtain the modified ion diameters alin
Coul , which
differ significantly from the ionic Born diameters aB . This demonstrates that the ionic Born diameters aB , which result from fitting thetotal solvation free energy to the Born equation, as given in Eq. (5),are not physically meaningful since the Born solvation model is notstrictly valid even for monovalent ions. The large difference between
the ionic Born diameters aB and the interaction ionic diameters d

In Fig. 11, we show the comparison between MSA theory
and GCMC simulations for dielectric decrements corresponding to
NaCl and NaI for various effective ionic diameters d. The agreement between theory and simulations for NaI using the diameterd = 0.27 nm (which is the best-fit result from Fig. 3) is not perfectfor concentrations above 2 M. The validity of MSA theory is governed by the electrostatic coupling parameter, which is the Coulomb
energy at contact, given by e2 /(4πε(cb )ε0 kB Td). For a concentrationof 2 M, the dielectric constant of NaI is rather small and given byε(2M) = 52.7, which results in a coupling parameter of 3.93. We
thus conclude that for coupling parameters below four, MSA theory reproduces well the GCMC simulation results. For NaCl and thebest-fit diameter d = 0.249 nm, the deviations between MSA and
GCMC are less because the dielectric constant is higher, and thus,
the electrostatic coupling parameter is smaller.
APPENDIX D: PREDICTIONS FOR ACTIVITY
COEFFICIENTS NEGLECTING THE IONIC SOLVATIONCONTRIBUTION
An interesting and somewhat puzzling observation is that theprimitive model reproduces the experimental activity coefficient
data very well if one neglects the ionic solvation contribution βμiw

FIG. 11. Mean ionic activity coefficient
as a function of salt concentration for (a)NaCl and (b) NaI. The lines represent
MSA theory combined with the Bornsolvation model [Eq. (2)] including thesalt-concentration-dependent dielectric
function. Symbols represent GCMC
simulation results using the same
parameters.

J. Chem. Phys. 153, 034103 (2020); doi: 10.1063/5.0016103
Published under license by AIP Publishing

153, 034103-12

The Journalof Chemical Physics

ARTICLE

scitation.org/journal/jcp

FIG. 12. Mean ionic activity coefficient as a function of salt concentration for (a) NaF, NaCl, and NaI, and (b) KF, KCl, and KI. The experimental data are represented by fullsymbols.83 Here, the theoretical predictions neglect the ionic solvation contribution βμiw in Eq. (4) to the mean ionic activity coefficient in the MSA result [Eq. (2)] and also in
the GCMC simulations. Empty symbols represent the GCMC simulations including fitted ionic diameters d and dielectric molten-salt parameters εms . The solid lines represent
MSA theory [Eq. (2)] including the salt-concentration-dependent dielectric function, while the dashed lines represent MSA theory for the same ion diameters d but for a fixeddielectric constant corresponding to pure water. The units were transformed from molality to molarity using experimental solution mass density data for KCl84 and NaCl.85 As
the correction is negligible for KCl and NaCl, we assume the same to be true also for the other salts (for which no data are available).

in Eq. (4) to the mean ionic activity coefficient in Eq. (2). In Fig. 12,
we compare the experimental data with MSA theory predictions aswell as GCMC simulations, where as the only fitting parameter wedetermine the effective mean ionic diameter d. Figure 12 is analogous to Fig. 3 where we do include the ionic solvation contribution
βμiw in both the MSA theory and the GCMC simulations. The fitted
ionic diameters d in Fig. 12 are different and significantly larger thanin Fig. 3. Curiously, the agreement between the experimental data
and the MSA and GCMC predictions is much better. In addition,
the agreement between the MSA and GCMC predictions is better,
which is due to the fact that the fitted diameter values in Fig. 12 arehigher.

DATA AVAILABILITY
The data that support the findings of this study are available
within the article.

REFERENCES
1

A. Heydweiller, Ann. Phys. 338, 145 (1910).
G. Jones and W. A. Ray, J. Am. Chem. Soc. 59, 187 (1937).
3
N. L. Jarvis and M. A. Scheiman, J. Phys. Chem. 72, 74 (1968).
4
P. K. Weissenborn and R. J. Pugh, J. Colloid Interface Sci. 184, 550 (1996).
5N. Matubayasi, K. Tsunetomo, I. Sato, R. Akizuki, T. Morishita, A. Matuzawa,
and Y. Natsukari, J. Colloid Interface Sci. 243, 444 (2001).
6
N. Matubayasi, K. Yamamoto, S.-i. Yamaguchi, H. Matsuo, and N. Ikeda,
J. Colloid Interface Sci. 214, 101–105 (1999).
7
W. Kunz, L. Belloni, O. Bernard, and B. W. Ninham, J. Phys. Chem. B 108, 2398
(2004).

2

J. Chem. Phys. 153, 034103 (2020); doi: 10.1063/5.0016103
Published under license by AIP Publishing

8

T. López-León, A. B. Jódar-Reyes, D. Bastos-González, and J. L. Ortega-Vinuesa,
J. Phys. Chem. B 107, 5696 (2003).
9
T. López-León, M. J. Santander-Ortega, J. L. Ortega-Vinuesa, and D. BastosGonzález, J. Phys. Chem. C 112, 16060 (2008).
10
J. M. Peula-García, J. L. Ortega-Vinuesa, and D. Bastos-González, J. Phys.
Chem. C 114, 11133 (2010).
11
A. P. dos Santos and Y. Levin, Phys. Rev. Lett. 106, 167801 (2011).
12
A. P. dos Santos, W. Figueiredo, and Y. Levin, Langmuir 30, 4593 (2014).
13
N. Jiang, P. Li, Y. Wang, J. Wang, H. Yan, and R. K. Thomas, J. Phys. Chem. B
108, 15385 (2004).
14
L. S. Romsted, Langmuir 23, 414 (2007).
15
J. Luczak, M. Markiewicz, J. Thöming, J. Hupka, and C. Jungnickel, J. Colloid
Interface Sci. 362, 415 (2011).
16
W. Müller, C. Déjugnat, T. Zemb, J.-F. Dufrêche, and O. Diat, J. Phys. Chem. B
117, 1345 (2013).
17
J. Morag, M. Dishon, and U. Sivan, Langmuir 29, 6317 (2013).
18
T. Oncsik, G. Trefalt, M. Borkovec, and I. Szilagyi, Langmuir 31, 3799 (2015).
19
Y. von Hansen, I. Kalcher, and J. Dzubiella, J. Phys. Chem. B 114, 13815 (2010).
20
H. I. Okur, J. Hladílková, K. B. Rembert, Y. Cho, J. Heyda, J. Dzubiella, P. S.
Cremer, and P. Jungwirth, J. Phys. Chem. B 121, 1997 (2017).
21
J. Lyklema, Chem. Phys. Lett. 467, 217 (2009).
22
W. Kunz, Curr. Opin. Colloid Interface Sci. 15, 34 (2010).
23
Y. Zhang and P. S. Cremer, Annu. Rev. Phys. Chem. 61, 63 (2010).
24
G. Trefalt, T. Palberg, and M. Borkovec, Curr. Opin. Colloid Interface Sci. 27, 9(2017).
25
M. Boström, D. R. M. Williams, and B. W. Ninham, Langmuir 17, 4475 (2001).
26
L. Perera and M. L. Berkowitz, J. Chem. Phys. 95, 1954 (1991).
27
L. X. Dang and D. E. Smith, J. Chem. Phys. 99, 6950 (1993).
28
P. Jungwirth and D. J. Tobias, J. Phys. Chem. B 105, 10468 (2001).
29
P. Jungwirth and D. J. Tobias, J. Phys. Chem. B 106, 6361 (2002).
30
P. Jungwirth and D. J. Tobias, Chem. Rev. 106, 1259 (2006).
31
D. Horinek and R. R. Netz, Phys. Rev. Lett. 99, 226104 (2007).
32
Y. Levin, Phys. Rev. Lett. 102, 147803 (2009).

153, 034103-13

The Journalof Chemical Physics

33

Y. Levin, A. P. dos Santos, and A. Diehl, Phys. Rev. Lett. 103, 257802 (2009).
D. J. Tobias, A. C. Stern, M. D. Baer, Y. Levin, and C. J. Mundy, Annu. Rev.
Phys. Chem. 64, 339 (2013).
35
Y. Levin and A. P. dos Santos, J. Phys.: Condens. Matter 26, 203101 (2014).
36
D. Horinek, A. Herz, L. Vrbka, F. Sedlmeier, S. I. Mamatkulov, and R. R. Netz,
Chem. Phys. Lett. 479, 173 (2009).
37
N. Schwierz, D. Horinek, and R. R. Netz, Langmuir 26, 7370 (2010).
38
S. I. Mamatkulov, C. Allolio, R. R. Netz, and D. J. Bonthuis, Angew. Chem., Int.
Ed. 56, 15846 (2017).
39
K. S. Pitzer, J. Phys. Chem. 77, 268 (1973).
40
K. S. Pitzer and G. Mayorga, J. Phys. Chem. 77, 2300–2308 (1973).
41
E. Moggia and B. Bianco, J. Phys. Chem. B 111, 3183 (2007).
42
C. Held, L. F. Cameretti, and G. Sadowski, Fluid Phase Equilib. 270, 87 (2008).
43
D. Fraenkel, J. Chem. Theory Comput. 11, 178 (2015).
44
L. Blum, Mol. Phys. 30, 1529 (1975).
45
H. R. Corti, J. Phys. Chem. 91, 686 (1987).
46
T. S. Sorensen, J. B. Jensen, and P. Sloth, J. Chem. Soc., Faraday Trans. 1 85,
2649 (1989).
47
J. B. Hasted, D. M. Ritson, and C. H. Collie, J. Chem. Phys. 16, 1 (1948).
48
W. R. Fawcett and A. C. Tikanen, J. Phys. Chem. 100, 4251 (1996).
49
J.-P. Simonin, L. Blum, and P. Turq, J. Phys. Chem. 100, 7704 (1996).
50
J.-P. Simonin, J. Phys. Chem. B 101, 4313 (1997).
51
I. Kalcher and J. Dzubiella, J. Chem. Phys. 130, 134507 (2009).
52
J. Vincze, M. Valiskó, and D. Boda, J. Chem. Phys. 133, 154507 (2010).
53
M. Valiskó and D. Boda, J. Chem. Phys. 140, 234508 (2014).
54
A. Levy, M. Bazant, and A. Kornyshev, Chem. Phys. Lett. 738, 136915 (2020).
55
A. Levy and D. Andelman, Phys. Rev. Lett. 108, 227801 (2012).
56
N. Gavish and K. Promislow, Phys. Rev. E 94, 012611 (2016).
57
C.-H. Ho, H.-K. Tsao, and Y.-J. Sheng, J. Chem. Phys. 119, 2369 (2003).
58
D. Henderson, D. Gillespie, T. Nagy, and D. Boda, Mol. Phys. 103, 2851
(2005).
59
A. Diehl, A. P. dos Santos, and Y. Levin, J. Phys.: Condens. Matter 24, 284115(2012).
60
L. Onsager and N. N. T. Samaras, J. Chem. Phys. 2, 528 (1934).
61
R. H. Garrett and C. M. Grisham, Biochemistry (Brooks/Cole, Boston, 2010).
62
Y. Marcus, J. Chem. Soc., Faraday Trans. 87, 2995 (1991).
63
P. Loche, C. Ayaz, A. Schlaich, D. J. Bonthuis, and R. R. Netz, J. Phys. Chem.
Lett. 9, 6463 (2018).
64
A. Bakhshandeh, A. P. dos Santos, A. Diehl, and Y. Levin, J. Chem. Phys. 142,
194707 (2015).
34

J. Chem. Phys. 153, 034103 (2020); doi: 10.1063/5.0016103
Published under license by AIP Publishing

ARTICLE

scitation.org/journal/jcp

65

M. P. Allen and D. J. Tildesley, Computer Simulations of Liquids (Oxford
University Press, Oxford, 1987).
66
Y. Uematsu, D. J. Bonthuis, and R. R. Netz, J. Phys. Chem. Lett. 9, 189 (2018).
67
D. Boda, D. Gillespie, W. Nonner, D. Henderson, and B. Eisenberg, Phys. Rev. E
69, 046702 (2004).
68
S. Tyagi, M. Süzen, M. Sega, M. Barbosa, S. S. Kantorovich, and C. Holm,
J. Chem. Phys. 132, 154112 (2010).
69
V. Jadhao, F. Solis, and M. de la Cruz, Phys. Rev. Lett. 109, 223905 (2012).
70
J. W. Zwanikken and M. Olvera de la Cruz, Proc. Natl. Acad. Sci. U. S. A. 110,
5301 (2013).
71
Z. Gan, H. Wu, K. Barros, Z. Xu, and E. Luijten, J. Comput. Phys. 291, 317(2015).
72
A. P. dos Santos, M. Girotto, and Y. Levin, J. Chem. Phys. 147, 184105(2017).
73
A. P. dos Santos and Y. Levin, J. Chem. Phys. 142, 194104 (2015).
74
A. P. dos Santos, M. Girotto, and Y. Levin, J. Phys. Chem. B 120, 10387 (2016).
75
P. W. Debye and E. Hückel, Phys. Z 24, 185 (1923).
76
Y. Levin, Rep. Prog. Phys. 65, 1577 (2002).
77
R. Buchner, G. T. Hefter, and J. Barthel, J. Chem. Soc., Faraday Trans. 90, 2475
(1994).
78
R. Buchner, G. T. Hefter, and P. M. Ma, J. Phys. Chem. A 103, 1 (1999).
79
J. Barthel, J. Krüger, and E. Schollmeyer, Z. Phys. Chem. 104, 59 (1977).
80
A. V. Kobelev, A. S. Lileev, and A. K. Lyashchenko, Russ. J. Inorg. Chem. 56,
652 (2011).
81
E. Glueckauf, Trans. Faraday Soc. 60, 1637 (1964).
82
A. Chandra, J. Chem. Phys. 113, 903 (2000).
83
W. J. Hamer and Y. C. Wu, J. Phys. Chem. Ref. Data 1, 1047 (1972).
84
M. T. Zafarani-Moattar and F. Izadi, J. Chem. Thermodyn. 43, 552 (2011).
85
P. S. Z. Rogers and K. S. Pitzer, J. Phys. Chem. Ref. Data 11, 15 (1982).
86
V. Balos, M. Bonn, and J. Hunger, Phys. Chem. Chem. Phys. 19, 9724 (2017).
87
N. Schwierz, D. Horinek, U. Sivan, and R. R. Netz, Curr. Opin. Colloid Interface
Sci. 23, 10 (2016).
88
J. M. Barthel, H. Krienke, and W. Kunz, Physical Chemistry of Electrolyte
Solutions: Modern Aspects (Springer, New York, 1998).
89
K. Ali, A. A. Shah, S. Bilal, and A. A. Shah, Colloids Surf., A 337, 194 (2009).
90
L. Piatkowski, Z. Zhang, E. H. Backus, H. J. Bakker, and M. Bonn, Nat.
Commun. 5, 4083 (2014).
91
T. Markovich, D. Andelman, and R. Podgornik, J. Chem. Phys. 142, 044702(2015).
92
D. E. Smith and L. X. Dang, J. Chem. Phys. 100, 3757 (1994).

153, 034103-14