Giant Axial Dielectric Response in
Water-Filled Nanotubes and Effective
Electrostatic Ion-Ion Interactions from a
Tensorial Dielectric Model
Philip Loche,† Cihan Ayaz,† Alexander Schlaich,‡ Yuki Uematsu,¶,§ andRoland R. Netz∗,††Fachbereich Physik, Freie Universität Berlin, 14195 Berlin, Germany‡Université Grenoble Alpes, CNRS, LIPhy, 38000 Grenoble, France¶Laboratoire de Physique de l’Ecole normale supérieure, ENS, Université PSL, CNRS,
Sorbonne Université, Université de Paris, F-75005 Paris, France§Department of Physics, Kyushu University, 819-0395 Fukuoka, Japan
E-mail: rnetz@physik.fu-berlin.de

Abstract

Introduction

Molecular dynamics simulations in conjunction
with effective medium theory are used to investigate dielectric effects in water-filled nanotubes. The resulting effective axial dielectric
constant shows a divergent increase for small
nanotube radii that depends on the nanotube
length, while the effective radial dielectric constant decreases significantly for thin nanotubes.
By solving Poisson’s equation for an anisotropic
dielectric medium in cylindrical geometry, weshow that the axial ion-ion interaction depends
for small separations primarily on the radial dielectric constant, not on the axial one. Thismeans that electrostatic ion-ion interactions inthin water-filled nanotubes are on the linear
dielectric level significantly enhanced due towater confinement effects at small separations,while at large separations the outside medium
dominates. If the outside medium is metallic,the ion-ion interaction decays exponentially forlarge ion separation.

Aqueous pores have key functional roles in cells,
membranes and proteins 1,2 , and are under intense investigation, for example in ion channels
and photosynthesic complexes 3,4 . Additionally,
they are interesting for technical applications
such as water desalination 5,6 , signal multiplication 7 and supercapacitors 8–10 A simple model
system for all these scenarios are water-filledcarbon nanotubes (CNTs), which are studied
intensively, both experimentally and in simulations 11–18 . In many applications, solvatedions are present in nanotubes and determine
most of the salient CNT properties. Consequently, ionic transport 19–21 and ionic conductivity 20,22,23 of CNTs are at the focus of investigations. To describe the static and kinetic ion
properties inside CNTs, the continuum PoissonBoltzmann and Helmholtz-Smoluchowski models are used 24–26 , but these models need as precise input predictions for the electrostatic ionion interactions, which in turn depend on the
dielectric properties. The water dielectric response close to planar interfaces has been shown
to be anisotropic and inhomogeneous 27–29 . In
particular, the perpendicular water dielectric1constant between two surfaces decreases significantly as the surface separation reaches thenanometer scale 30,31 . Also for water confined
inside CNTs the dielectric properties have beenshown to differ from bulk water using simulations 32,33 as well as experiments 34,35 . However,to derive from dielectric properties the electrostatic interaction between ions requires the solution of the Poisson equation in the presenceof tensorial and inhomogeneous dielectric profiles in cylindrical symmetry, which is far fromtrivial.
From molecular dynamics simulations of water inside single-walled armchair carbon nanotubes of various radii and lengths, we extract
the tensorial dielectric radial and axial profiles. The definition of these profiles reflects thelong range and non-locality of the dielectric response 36,37 and has been intensely discussed inliterature 38,39 . With decreasing CNT radius, wefind the axial dielectric response to be dramatically enhanced and to become CNT length dependent, which reflects the long-ranged dipolarordering of water. The inverse radial dielectric
response shows oscillations and multiple signchanges. For spatially varying dielectric tensorial profiles the electrostatic Green’s function,i.e., the interaction between two point charges,
cannot be calculated in closed form. We, therefore, use effective medium theory to extracteffective axial and radial dielectric constants.
These dielectric constants become anisotropic
for nanotube radii R < 10 , nm, with the axial value diverging in a CNT-length dependent
fashion and the radial value approaching thevacuum value in the single-file limit. By comparison with corresponding results for water inplanar confinement, we show that the confinement effects in cylindrical confinement are evenmore drastic than in planar confinement, whichsuggests that curvature plays an important rolein confinement. We also establish the Green’s
function formalism for tensorial electrostatics
in cylindrical confinement. We find that for
CNTs the axial ion-ion interaction is governedby the outside dielectric constant for large ionion distances, we present explicit results for thecases where the outside medium is dielectricallyequivalent to vacuum and where the outside

medium is metallic. For very small separations,the interaction only depends on the radial dielectric constant, not on the axial one, and thus
is dramatically enhanced compared to the ionion interaction in bulk water. The axial ion-ion
interaction in water-filled CNTs thus depends
primarily on the radial dielectric water properties and exhibits complex crossovers. The most
relevant message for experiments with ions in
CNTs is the crossover due to the reduced radialdielectric component for small separation andthe influence of the outside dielectric constantfor large separations, which will influence ioniccorrelations and transport properties as well as
chemical reactions 40 .

Methods
TheoryFor a cylindrical system that is homogenous
along the axial coordinate z and the angularcoordinate ϕ, the ϕ component of the rotationof the static electric field follows from Maxwell’sequations as∂Er (r) ∂Ez (r)
−
= 0.
∂z
∂r

(1)

Since the electric field is independent of z, thefirst term vanishes, and thus Ez (r) = const. If
no free charges are present, Gauss’s law for the
displacement field reads
1 ∂rDr (r) = 0
r ∂r

(2)

and thus rDr (r) = const. In linear response
theory, the change of the displacement field,
D(r), is related linearly to the change of theelectric field, E(r), via the non local dielectric
tensor ε̂nl (r, r 0 ) according toZ
D(r) = ε0ε̂nl (r, r 0 )E(r 0 )dr 0 ,(3)Vwhere V is the system volume. Using the condition that Ez (r) is constant, as follows from eq

2

1, the z component of eq 3 simplifies toternal field. In this expression dX denotes thephase space integral over all positional degrees
of freedom. Using eq 10 in eq 8, we find for theaxial dielectric profile


Z ∞
β2πL
0 0
0
mz (r )r drmz (r)εz (r) = 1 +
ε0
0

Z ∞
0 0
0
.
− hmz (r)imz (r )r dr

Dz (r) = ε0 εz (r)Ez
(4)
R
0
0
with εz (r) = V εzz
Similarly, thenl (r, r )dr .
electric field can be written in terms of thechange of the displacement field and the inverse
dielectric response function as
Z
−1
0
0
0
E(r) = ε0
ε̂−1(5)nl (r, r )D(r )dr .

0

(11)

V

Similarly, by using eq 10 in eq 9 and employing
eq 7, we obtain for the radial dielectric profile


Z ∞
β2πrL
−1
0
0
εr (r) = 1 −mr (r )drmr (r)
ε0
0
Z ∞

0
0
− hmr (r)i
mr (r )dr.

Using the condition that rDr (r) is constant, as
follows from eq 2, the r component of eq 5 isfound to be−1
Er (r) = ε−10 εr (r)Dr (r) ,(6)

R ∞ rr0 −1where ε−1
dr0 . Byr (r) = 2πrL 0 [εnl (r, r )]inserting the relation
D(r) = ε0 E(r) + m(r)

0

(12)

Note that no locality assumptions are used toderive eqs 11 and 12, which agree with previous results obtained by slightly different arguments 33 .

(7)

into eq 4, where m(r) denotes the difference ofthe polarization density with and without anapplied external field, one obtains for the linear
response axial dielectric profile
εz (r) = 1 +Computational Details

1 ∂mz (r)mz (r)≈1+,
ε0 Ez
ε0 ∂Ez Ez =0
(8)

(a)

(b)

where the response to the electric field is included to linear order by expanding mz (r). The
linear response radial dielectric profile is obtained by substituting eq 7 into eq 6
ε−1
r (r) = 1 −∂mr (r)mr (r)≈1−r.
Dr (r)∂(rDr (r)) rDr (r)=0(9)

Figure 1: (a) Simulation snapshot of a waterfilled CNT with radius R = 1.2 nm. (b) Representation of the effective medium model used
To obtain an expression for the polarization
for solving Poisson’s equation. The permittivdensity derivative in eqs 8 and 9, we use theity is constant and anisotropic inside the Gibbs
canonical average
dividing radius RGDS and isotropic and unity



R
R
outside.
m(r) exp −β H − E(r 0 )m(r 0 )dr 0 dX



R
R
hm(r)i =,exp −β H − E(r 0 )m(r 0 )dr 0 dX We simulate a periodic single-walled armchair(10)
CNT filled with water in a rectangular simulation box with periodic boundary conditions.
where β is the inverse thermal energy and H is
The carbon atoms are frozen, unpolarizable andthe system Hamiltonian in the absence of an excharge-neutral, for water we use the SPC/E
3

The axial permittivity profiles εz (r) in Figure
2b are similar in shape to the density profiles,
but for smaller radii R the maxima in the εz (r)profile become increasingly pronounced. The
inverse radial dielectric profiles ε−1
r (r) in Figure 2c are very different from εz (r) and exhibit
pronounced oscillations in the CNT center. To
clarify the origin of these oscillations we extractε−1
r (r) from simulations of a large homogeneous
water box with periodic boundary conditions.
The result, shown as a black dashed line in Figure 2c, perfectly superimposes on the data for
R = 1.7 nm and R = 2.7 nm. This shows thatthese oscillations are not caused by the presence of the CNT but rather are due to the water discreteness. In Figure 2d–f, the profiles are
shifted by the position of the nanotube radius
R. It becomes clear that the oscillations in n(r)and εz (r) (and to a lesser degree also in ε−1r (r))
close to the CNT interface are rather independent of R and agree quite nicely with the planar
limit (solid black line), obtained from simulations of water between planar graphene sheetswith a large separation. Only the axial dielectric profile for the narrowest CNT (R = 0.4 nm)
is length dependent, as shown in Figure S4 anddiscussed below.

model, see Figure 1a for a simulation snapshot.
We define the CNT radius R by the carbon
atom position as given in eq S1 in the Supporting Information and shown in Figure 1b.
For the carbon-oxygen Lennard-Jones interaction, we use the GROMOS53a6 force field 41 , asit accurately reproduces the water contact angle on carbon surfaces 42 . The dielectric profileat planar surfaces has been previously shown todepend slightly on the surface-water interaction
strength 28 , we expect a similar dependence also
for the present cylindrical case, which however
is not studied further. To determine the correct
water density inside the CNTs, we first attacha reservoir to cylindrical tubes with a lengthof L = 15.1 nm, from which we determine thewater density in the central CNT section. The
reservoir is pressure-coupled along the cylinder
axes using an extended Hamiltonian and the
Berendsen method, see section S2 of the Supporting Information for details. For the production runs, we use periodic CNTs withouta reservoir in the NVT ensemble and the vrescale thermostat 43 at 300 K. Different temperature control methods have no effect on our
results as shown in section S3 of the Supporting
Information. Simulations were performed using the GROMACS simulation package 44 andwe provide a free Python library (located athttps://gitlab.com/netzlab/maicos) basedon the MDAnalysis library 45 , which implements
the analysis we used for extracting the dielectric
profiles 28,30,39 . See sections S1 in the Supporting Information for further simulation details.

Effective Medium Theory
In order to calculate ion-ion electrostatic interactions in CNTs analytically, the r-dependent
dielectric tensorial profiles must be modeledby piecewise constant profiles using effective
medium theory 46 , as schematically indicated in
Figure 1b. This construction is subtle and requires care, to make this mapping unique werequire asymptotic integral invariances that are
different for the axial and radial components.
For the axial part we demand that the laterally averaged displacement field Dz in the
CNT is the same for the actual dielectric profile εz (r) and the effective piece-wise constant

Results and Discussion
MD Simulations
From the MD simulations, the radial number
density n(r) and permittivity profiles εz (r) and
ε−1r (r) are deduced according to eqs 11 and 12
and shown in Figure 2 for different CNT radii
R. Profiles are for CNT length L = 15.1 nm.
Different lengths are shown in Figure S4 andS5 in section S4 of the Supporting Information.
The density profiles n(r) in Figure 2a show thetypical oscillations close to the CNT wall 14 .

4

(b) 600

50

0

1

150

400
200
0

2

r / nm

(e)

(c)

0

600

1

0
−1.00 −0.75 −0.50 −0.25
r − R / nm

−20

0
(f)

400
200

0.00

0

2

r / nmε−1
r (r )εz ( r )n(r ) / nm−350periodic water box−40

gr aphene

100

20

0
−1.00 −0.75 −0.50 −0.25
r − R / nm

0.00

1

2

r / nm

2

50

1

25

0

0

−1

−25

−2
−1.00 −0.75 −0.50 −0.25
r − R / nm

ε−1
r (r )

R = 2.7 nm

ε−1
r (r )

1000(d)

R = 0.7 nm
R = 1.7 nm

εz ( r )n(r ) / nm−3(a) 150−50

0.00

Figure 2: (a) Water oxygen number density profiles n(r), (b) axial dielectric profiles εz (r), (c)inverse radial dielectric profiles ε−1
r (r) for CNTs of several radii and fixed length L = 15.1 nm. The
horizontal broken lines depict the bulk values obtained from independent simulations. (d–e) The
same results as in a–c, but shifted by the respective CNT radius R, together with the profile ata planar graphene interface, shown as a solid black line. The right axis in (f) applies to data for
R = 0.7 nm.
profile ε∗z (r). Using eq 4 we obtainrather independent of our choice for the radiusof the effective dielectric profile, as we show insection S5 in the Supporting Information.
For the radial component, we demand the potential difference between inside and outside the
CNT to be the same when using the actual dielectric profile εr (r) and the effective box profileε∗r (r). From the relation between the electrostatic potential Φ(r) and the electric field E(r),∂Φ(r)/∂r = −Er (r), using eq 6 and the factthat rDr (r) is constant, we obtain

Z R
1
2πrDz (r)dr
Dz =
πR2 0
Z
Z R
2Ez ε0 R! 2Ez ε0=rεz (r)dr =rε∗z (r)dr ,
2
R2
R
0
0
(13)
where the integral ranges up to the CNT radius
R. The effective profile ε∗z (r) is chosen to be abox profile with an effective dielectric constantof εeffz ,(εeff
if r ≤ RGDSzε∗z (r) =.
(14)
1
if r > RGDS

Φ(R) − Φ(R0 ) = ε−1
0 rDr (r)!

= ε−10 rDr (r)Substituting this into eq 13, yields εeff
z as
εeff
z −1 =

RR

2 0 r (εz (r) − 1) dr(RGDS )2

.

Z R

1

0
0
R0 r εr (r )
Z R

1

0 ∗ 0
R0 r εr (r )

dr0
dr0 .
(16)

As for the axial component, ε∗r (r) is a box profile(
1/εeff
if r ≤ RGDS
1
r
=
,
(17)
ε∗r (r)1
if r > RGDS

(15)

The
Gibbs dividing radius is RGDS =p(N vw )/(πL), where N is the number of water molecules and vw = 0.0307 nm3 is the bulk
molecular water volume. The result for εeff
z iswhere εeff
r is the effective radial dielectric con5

stant. Substituting eq 17 into eq 16, εeff
r is determined asRR(εr (r)r)−1 dr − ln (R/R0 )1
R0. (18)=εeff
ln (RGDS /R0 )rdering and agrees nicely with analytical models
for single-file water chains 47 , as we show in Figure S6 in the Supporting Information.
This anisotropic behavior of the dielectric behavior is somewhat similar to results for planar
confinement 30 , where the dielectric component
perpendicular to the surface was found to decrease and the parallel component was foundto increase slightly with decreasing surface separation. To compare dielectric confinement effects in cylindrical to planar systems we calculate effective dielectric constants of water between two planar graphene sheets at a separation L = 2R using the scheme introduced inRef. 30. From the data shown in Figure 3, weconclude that the dielectric behavior of waterin CNTs is rather different compared to water between planar graphene sheets and that
confinement effects in cylinders are more pronounced than between planar sheets. This suggests that curvature effects play an importantrole in confinement. In fact, similar curvatureeffects were previously found also for the water
friction coefficient in CNTs 17 . It is well known
that the dielectric constant in bulk decreases
significantly with rising temperature 48 , an effect that is well reproduced in simulations 49 .
In section S7 of the Supporting Information weshow that the axial component of the effective
dielectric constant decreases with temperature
whereas the radial effective dielectric constantincreases with temperature. Thus, the temperature dependence in cylindrical confinement
differs from the bulk behavior, which is an interesting finding.

The lower integration boundary R0 in eq 18
corresponds to the radius at which the potential
inside the CNT is determined. It is chosen suchthat εeff
r for the largest radius R = 10.17 nm
matches the SPC/E bulk permittivity of ε = 70.
There are several choices for R0 that comply
with this condition, we pick R0 = 0.38 nm forwhich εeff
r grows monotonically with R (details
are given in section S6 in the Supporting Information).
105
L = 5.1 nm
L = 15.1 nm
L = 75.1 nm
L = 150.1 nm

104

εeff

103

εeff
εeff
r
z

102
εeff
⊥ , graphene L = 2R

101

εeff
k , graphene L = 2R

100
0

2

4

6
R / nm8

10

Figure 3: Effective radial (circles) and axial
and(triangles) dielectric water constants εeffr
eff
εz for several CNT lengths and radii. The horizontal dotted line shows the bulk dielectric constant of SPC/E water ε = 70. Open black circles and triangles show the effective perpendicular and parallel dielectric constants between
planar graphene sheets at separation of L = 2R.

Collective Versus Self Polarization

Figure 3 shows the effective dielectric coneff
stants εeffz and εr calculated according to eqs
15 and 18. The radial part decreases whereas
the axial part diverges as R → 0. This behavior
reflects the dipolar ordering of water moleculesalong the CNT axis, which becomes stronger
with decreasing CNT radius, as will be explained by a simple Langevin model further below. For the smallest CNT radius, R = 0.4 nm,
εeff
z in fact increases dramatically with increasing CNT length. This is caused by the increasing correlation length of the dipolar water or-

In order to understand the cause of the pronounced difference between axial and radial dielectric effects, we investigate self and collective
polarization fluctuations. For this we define theself and collective polarization variances asX
Cαself =hp2i,α i − hpi,α i2(19)iCαcoll =

XX
i

6

j6=i

hpi,α pj,α i − hpi,α ihpj,α i ,(20)where pi,α is the dipole component α ofmolecule i and α = x, y, z.
We definethe average radial components as Crself =(Cxself + Cyself )/2 and Crcoll = (Cxcoll + Cycoll )/2.
The total polarization variance Cα is given by
Cα = Cαself + Cαcoll , and in bulk the relation
between Cα and the dielectric constant reads 37
εα = 1 +

103
model
model

Cα / ε0 π R GDS

2

LkB T

(a)

102

CzselfCzcollCrselfCrcoll

101

g

gz

101

gr

100
10−1

0

1

2

3
R / nm4(21)

We note that there is a subtle difference between Cr and the effective radial dielectric constant εeff
r , as an inspection of the defining eq
12 for the radial dielectric profile shows. Figure 4a shows the polarization variances rescaled
GDS 2
by ε0 π R
LkB T , so that the data can becompared to the dielectric constants calculated
earlier. For the axial part, collective contributions are positive and much larger than theself part. As shown in Figure 4b, the axial
Kirkwood gz factor, defined as gz = Cz /Czself ,reaches for large CNT radii a value of gz ≈3, close to what is found in bulk water 50 ,but for smaller radii drastically increases. For
larger CNT radii therefore roughly three water molecules are axially correlated, while thisnumber rapidly increases as the radius goes
down. For the radial contribution, the behavior is completely different. Here, the collective
contribution is negative, which means that theradial component of the orientation of adjacentwaters is anti-correlated. We note that for verylarge radii the radial and axial collective contributions must become equal, but this crossoveris very slow. In accordance with this, the radial
Kirkwood factor gr = Cr /Crself in Figure 4b isvery small and approaches the bulk value very
slowly.

−101
(b)

Cα
.
ε0 V kB T

5

Figure 4: (a) Self and collective contributions of
the rescaled polarization variances for axial (triangles) and radial (circles) directions. Solid andbroken lines denote results from the anisotropic
Langevin model eqs 24 and 25 for the self polarization variances. The horizontal dotted line
shows the bulk dielectric constant of SPC/E
water ε − 1 = 69. (b) Kirkwood factors gz andgr obtained from axial and radial polarization
fluctuations. The horizontal dotted line shows
the bulk water Kirkwood factor g = 2.68 50 .

Anisotropic Langevin Model
Here we generalize the Langevin model forthe response of a single dipole to an applied
electric field to account for anisotropic alignment effects. The Hamiltonian for the standard
Langevin model reads H = −p · E. We extend
this by an biaxial alignment field with strength7Electrostatic Interactions in an
Anisotropic Cylindrical Dielectric
H = −p [Ez cos(ϑ) + Ex sin(ϑ) cos(ϕ)] + c cos(ϑ) ,Medium(22)
c according to

According to the effective medium model we derived from our simulation data, inside a cylinder of radius RGDS the dielectric constants areeffanisotropic and given by εeffz and εr , while outside the cylinder the dielectric constants are different and reflect the dielectric properties of theembedding medium, see Figure 1b. The electrostatic potential Φ(r) for an arbitrary charge
distribution ρ(r) follows from Poisson’s equation

which mimics the axial orientation of the water dipole along z in narrow CNTs for c < 0.
We assume an external electric field in thex − z plane and describe the dipole orientation by spherical coordinates with polar and azimuthal angles ϑ and ϕ. The permittivities are
obtained from second derivatives of the partition function
Z, the integral over phase spaceR
Z = dX exp(−βH), with respect to the applied field E



∂
1 ∂Z
n
, (23)
εα = 1 +βε0 ∂Eα Z ∂Eα Eα =0ε0 {∇ · [ε̂(r) · ∇Φ(r)]} = −ρ(r) ,(26)

where ε̂(r) is the anisotropic permittivity tensor. The derivation of Φ(r) for a point charge
distribution ρ(r) = qδ(r − r 0 ) is explicitly
shown in section S9 in the Supporting Information.
We consider the interaction energy U (z) =q 0 Φ(z) between two point charges that are located at the CNT axis at separation z. Fromthe exact expression for U (z), eq S50 in the
Supporting Information, we find for small distances z the asymptotic behavior


1 α−1
γR
qq 0
+
,
lim U (z) =z→0
4πε0 εeffz α + 1 (z/2)2 + (γR)2r(27)

where α = z, x. For the two components ofthe permittivity, as shown in section S8 of the
Supporting Information, we find
√


np2
2 βce−βc√εz = 1 +
,
(24)
1− √
2cε0
πerf( βc)√
np2
2 βce−βc√εx = 1 +1 − 2βc − √,4cε0πerf( βc)(25)
where n is the number density and erf(x) theerror function. In the limit of c = 0 one recovers
the standard Langevin model result εz = εx =1 + (βnp2 )/(3ε0 ), which gives εz = εx = 19.71
for a dipole moment p = 0.049 e nm and n =33.4nm−3 for the SPC/E water model. Predictions for εz − 1 (from eq 24) and εx − 1 (fromeq 25) are included in Figure 4a as solid andbroken lines and describe the self-polarization
very well, where the radius-dependent alignment field strength c(R) = c0 + bR−n was fitted
to the simulation data as shown in Figure S13
in the Supporting Information. The best fit isobtained for c0 = −0.2 kB T, b = 0.28 kB T nmn
and n = 2.8. This shows that the symmetrybreaking between the axial and radial dielectricself responses can be described by an extended
Langevin model that accounts for a radiusdependent axial water alignment.

pp
eff and γ =
eff
ε
εeff
where α = εeff
z r
z /εr . For
small z the first term in eq 27 dominates andthe interaction energy only depends on the effective radial dielectric constant εeff
r but not on
the axial one. As we show in section S9.1 in
the Supporting Information, this at first sight
surprising behavior follows quite generally foranisotropic homogeneous dielectric media andis not specific to the cylindrical geometry. In
the limit z → ∞ we findlim U (z) =z→∞

qq 0 1.
4πε0 εout z(28)

As expected, for large z the interaction energy
only depends on the outer dielectric constantεout . For convenient use in coarse-grained sim8

ulations we provide a simple heuristic crossover
formula for U (z) that exhibits the correct
asymptotic behaviors in section S9.5 in the Supporting Information.

U (z ) / qq 0 kB T(a)

U (z ) / qq 0 kB T(b)

ent CNT radii as solid colored lines. For thedashed colored lines in Figure 5a we also use
eq S50 but rather take isotropic dielectric constants that correspond to SPC/E bulk waterεeff= εeff= 70 inside the CNTs. The difrzference between these two predictions indicatesthat the interaction is in CNTs considerably enhanced below separations z < 2 nm due to thedecrease of εeff
r compared to bulk water. We
also see that for distances z > 2 nm the interactions deviate from a pure power law U (z) ∼ 1/z
and thus are influenced by the outside medium.
For large distances z > 10 nm the interactions
only depend on the permittivity of the outside
medium and approach eq 28, which is shown asa solid black line. The dotted black line shows
the interaction in bulk water. We see that theinteractions inside the CNTs are for an intermediate distance range 10 nm< z < 100 nm
even larger than in vacuum, which is caused
by image charge effects. We conclude that theinteractions between charges in CNTs are enhanced compared to bulk water because of thepresence of a low-dielectric medium outside the
CNT and because the radial effective dielectric
constant is significantly reduced in small CNTs
compared to bulk water.
The case where the outside medium is metallic is of considerable interest since carbon nanotubes can be conducting 51 . In Figure 5b we
present results for the limiting case where theouter dielectric constant diverges, εout = ∞, bysolid lines. Due to the strong image charge effects, the interaction between charges becomes
exponentially screened for large z and decays as103
R = 0.7 nm
R = 1.7 nm
R = 2.7 nm

102
101
100
10−1
10−2
103

lim U (z )z →∞Ubulk (z )

102
101
100
10−1
10−2

10−3

lim U (z )

z →∞

10−2

10−1

100
z / nm

101

102

103

Figure 5: Electrostatic interaction energy U (z)between two charges as a function of the distance z along the axial direction inside a cylinder that is filled with an anisotropic dielectric
medium, the outside medium in (a) is an insulator and modeled as vacuum and in (b) is aconductor and modeled as a metal. (a) Coloredlines show U (z) according to the exact solutionof the Poisson equation, eq S50 in the Supporting Information. For the solid colored lines,
we use the effective anisotropic dielectric constants as shown in Figure 3. For the dashed colored lines we use isotropic dielectric constantseffof εeffr = εz = 70. The dotted black line shows
the ion-ion interaction in isotropic bulk waterwith ε = 70, the solid black line shows the ionion interaction in vacuum with ε = 1. (b) Solidlines show the exact solution inside a metallic
medium according to eq S71 in the Supporting
Information. The dashed dotted lines are thelimiting result for large z according to eq. 29.

lim U (z) =z→∞1z
qq 0 J2 (ξ1 ) − ξγRe,πε0 αR J2 (ξ1 )(29)

where Jn (z) is the n-th order Bessel function of
first kind and ξ1 ≈ 2.4048 is the first zero ofJ0 . The exact solution of the metallic limit andthe derivation of the limiting law is provided in
section S9.6 of the Supporting Information.
The horizontal gray lines in Figures 5 a and b
denote the Bjerrum length lB , which is the distance at which the electrostatic interaction between two unit charges equals kB T . The Bjerrum length inside CNTs with diameters of the

In Figure 5a we show results according to
the exact expression for U (z) given in eq S50
for the case where the outside medium is vacuum and thus εout = 0; here we use the effective anisotropic dielectric constants inside thecylinder presented in Figure 3 for three differ9

order of a few nanometers embedded in vacuum
is seen to be about 100 nm and thus two ordersof magnitude larger than in bulk water, where
we have lB ≈ 0.7 nm. For metallic CNTs the
Bjerrum length is about one order of magnitude larger compared to bulk water. This gives
the electrostatic interaction in CNTs an enormous range that will dramatically enhance correlation effects.

large separations and thus are enhanced compared to interactions in bulk water. If the outside medium is metallic the ion-ion interactions
become exponentially screened for large separations.
Clearly, for ions of finite size, we expect deviations from the present calculation for point-likeions. In fact, extended charge distributions canbe treated by Green’s function methods thatwere used earlier by us for planar systems 52 .
The effective anisotropic dielectric model forwater inside CNTs developed in this paper will
also be useful for future coarse-grained modeling of confined aqueous systems 53,54 . In particular, the long-range nature of electrostatic ionion interactions in water-filled cylinders will be
important to take into account in future modeling of ionic distributions, ionic conductivities
and electrokinetic effects inside CNTs.

Conclusions
We analyze the dielectric properties of waterfilled CNTs using MD simulations combinedwith effective medium theory. We find thatthe axial dielectric profile oscillates close to the
CNT wall with an amplitude that depends on
the CNT radius. This is due to water layeringand matches well with the oscillations in thewater density profiles. The behavior of the radial dielectric profiles is markedly different anddoes not reflect the water density profile. Effective medium theory is used to extract effective
dielectric constants from these dielectric profiles. This, on the one hand, helps the interpretation of the dielectric effects and, on the other
hand, enables us to calculate electrostatic interactions inside CNTs. We find that the effective
axial dielectric constant dramatically increases
with decreasing CNT radius, which by comparison with a simple extended Langevin model wetrace back to axial water ordering in thin CNTs.
The effective radial dielectric constant, in contrast, decreases with decreasing CNT radius.
By the exact solution of the anisotropic
Poisson equation in cylindrical geometry, wedemonstrate that the axial interaction between
two ions inside a CNT depends for small separation primarily on the radial dielectric constant,
not on the axial one. Thus, the ion-ion interaction is significantly enhanced for small separation due to the decrease of the radial dielectric constant for small CNT radii, compared tothe hypothetical scenario where the dielectric
constant would be isotropic and given by thebulk water value. If the outside medium has
a low dielectric constant, the ion-ion interactions are determined by the outside medium for

Acknowledgement
FundingbytheDeutsche Forschungsgemeinschaft (DFG) via
grants SFB 1078 and SFB 1114 is acknowledged.

Supporting Information Available
Simulation details, determination of the water
density inside carbon nanotubes, dielectric profiles for different CNT lengths, comparison of
the axial effective dielectric constant for singlefile water with an analytical model, thermostat
and temperature effects on the dielectric constant, the effective dielectric constant for different effective radial and axial radii, results fordifferent choices of the inner radius R0 , derivation of the extended Langevin model, derivationof the tensorial electrostatic model

References(1) Sui, H.; Han, B.-G.; Lee, J. K.; Walian, P.;
Jap, B. K. Structural basis of waterspecific transport through the AQP1 water channel. Nature 2001, 414, 872–878.
(2) Rasaiah, J. C.; Garde, S.; Hummer, G.
10

Water in Nonpolar Confinement: From
Nanotubes to Proteins and Beyond. Annu.
Rev. Phys. Chem. 2008, 59, 713–740.

(12) Majumder, M.; Chopra, N.; Andrews, R.;
Hinds, B. J. Nanoscale hydrodynamics:
Enhanced flow in carbon nanotubes. Nature 2005, 438, 44.

(3) Agre, P. The Aquaporin Water Channels.
Proc Am Thorac Soc 2006, 3, 5–13.

(13) Whitby, M.; Cagnon, L.; Thanou, M.;
Quirke, N. Enhanced Fluid Flow throughNanoscale Carbon Pipes. Nano Lett.
2008, 8, 2632–2637.

(4) Pfeil, B. E.; Schoefs, B.; Spetea, C. Function and evolution of channels and transporters in photosynthetic membranes.
Cell. Mol. Life Sci. 2014, 71, 979–998.

(14) Alexiadis, A.; Kassinos, S. Molecular Simulation of Water in Carbon Nanotubes.
Chem. Rev. 2008, 108, 5014–5034.

(5) Corry, B. Designing Carbon Nanotube
Membranes for Efficient Water Desalination. J. Phys. Chem. B 2008, 112, 1427–
1434.

(15) Hummer,
G.;
Rasaiah,J.
C.;
Noworyta, J. P. Water conduction
through the hydrophobic channel of acarbon nanotube. Nature 2001, 414,
188–190.

(6) Shannon, M. A.; Bohn, P. W.; Elimelech, M.; Georgiadis, J. G.; Mariñas, B. J.; Mayes, A. M. Science and technology for water purification in the coming
decades. Nature 2008, 452, 301–310.

(16) Mukherjee, B.; Maiti, P. K.; Dasgupta, C.;
Sood, A. K. Single-File Diffusion of Water Inside Narrow Carbon Nanorings. ACS
Nano 2010, 4, 985–991.

(7) Tu, Y.; Xiu, P.; Wan, R.; Hu, J.; Zhou, R.;
Fang, H. Water-mediated signal multiplication with Y-shaped carbon nanotubes.
PNAS 2009, 106, 18120–18124.

(17) Falk, K.; Sedlmeier, F.; Joly, L.;
Netz, R. R.; Bocquet, L. Molecular Origin
of Fast Water Transport in Carbon Nanotube Membranes: Superlubricity versus
Curvature Dependent Friction. Nano Lett.
2010, 10, 4067–4073.

(8) Frackowiak, E.; Jurewicz, K.; Delpeux, S.;
Béguin, F. Nanotubular materials for supercapacitors. Journal of Power Sources
2001, 97-98, 822–825.

(18) Chakraborty, S.; Kumar, H.; Dasgupta, C.; Maiti, P. K. Confined Water:
Structure, Dynamics, and Thermodynamics. Acc. Chem. Res. 2017, 50, 2139–2146.

(9) Raymundo-Piñero, E.; Kierzek, K.; Machnikowski, J.; Béguin, F. Relationship between the nanoporous texture of activatedcarbons and their capacitance properties
in different electrolytes. Carbon 2006, 44,
2498–2507.

(19) Samoylova, O. N.; Calixte, E. I.; Shuford, K. L. Molecular Dynamics Simulations of Ion Transport in Carbon Nanotube Channels. J. Phys. Chem. C 2015,
119, 1659–1666.

(10) Feng, G.; Li, S.; Presser, V.; Cummings, P. T. Molecular Insights into Carbon Supercapacitors Based on RoomTemperature Ionic Liquids. J. Phys.
Chem. Lett. 2013, 4, 3367–3376.

(20) Secchi, E.; Niguès, A.; Jubin, L.; Siria, A.;
Bocquet, L. Scaling Behavior for Ionic
Transport and its Fluctuations in Individual Carbon Nanotubes. Phys. Rev. Lett.
2016, 116, 154501.

(11) Holt, J. K.; Park, H. G.; Wang, Y.; Stadermann, M.; Artyukhin, A. B.; Grigoropoulos, C. P.; Noy, A.; Bakajin, O.
Fast Mass Transport Through Sub-2Nanometer Carbon Nanotubes. Science
2006, 312, 1034–1037.

(21) Secchi, E.; Marbach, S.; Niguès, A.;
Stein, D.; Siria, A.; Bocquet, L. Massive
radius-dependent flow slippage in carbon
nanotubes. Nature 2016, 537, 210–213.
11

(22) Zwolak, M.; Wilson, J.; Ventra, M. D. Dehydration and ionic conductance quantization in nanopores. J. Phys.: Condens.
Matter 2010, 22, 454126.

Watanabe, K. et al. Anomalously low dielectric constant of confined water. Science 2018, 360, 1339–1342.
(32) Mikami, F.; Matsuda, K.; Kataura, H.;
Maniwa, Y. Dielectric Properties of Water
inside Single-Walled Carbon Nanotubes.
ACS Nano 2009, 3, 1279–1287.

(23) Kavokine, N.; Marbach, S.; Siria, A.; Bocquet, L. Ionic Coulomb blockade as afractional Wien effect. Nat. Nanotechnol.
2019, 14, 573–578.

(33) Schaaf, C.; Gekle, S. Spatially resolved dielectric constant of confined water and its
connection to the non-local nature of bulk
water. The Journal of Chemical Physics
2016, 145, 084901.

(24) Schoch, R. B.; Han, J.; Renaud, P. Transport phenomena in nanofluidics. Rev.
Mod. Phys. 2008, 80, 839–883.
(25) Bocquet, L.; Charlaix, E. Nanofluidics,
from bulk to interfaces. 2010, 39, 1073–
1095.

(34) Cui, H.-B.; Takahashi, K.; Okano, Y.;
Kobayashi, H.; Wang, Z.; Kobayashi, A.
Dielectric Properties of Porous Molecular
Crystals That Contain Polar Molecules.
Angewandte Chemie International Edition
2005, 44, 6508–6512.

(26) Bonthuis, D. J.; Netz, R. R. Unravelingthe Combined Effects of Dielectric andViscosity Profiles on Surface Capacitance,
Electro-Osmotic Mobility, and ElectricSurface Conductivity. Langmuir 2012, 28,
16049–16059.

(35) Zhou, B.; Kobayashi, A.; Cui, H.-B.;
Long, L.-S.; Fujimori, H.; Kobayashi, H.
Anomalous Dielectric Behavior and Thermal Motion of Water Molecules Confinedin Channels of Porous Coordination Polymer Crystals. J. Am. Chem. Soc. 2011,
133, 5736–5739.

(27) Stern, H. A.; Feller, S. E. Calculationof the dielectric permittivity profile for anonuniform system: Application to a lipid
bilayer simulation. J. Chem. Phys. 2003,
118, 3401–3412.

(36) Fröhlich, H. Theory of Dielectrics: Dielectric Constant and Dielectric Loss; Oxford
University Press, 1987.

(28) Bonthuis, D. J.; Gekle, S.; Netz, R. R.
Profile of the static permittivity tensor ofwater at interfaces: consequences for capacitance, hydration interaction and ion
adsorption. Langmuir 2012, 28, 7679–
7694.

(37) Neumann, M. Dipole moment fluctuation
formulas in computer simulations of polar systems. Molecular Physics 1983, 50,
841–858.

(29) Zhang, C.; Gygi, F.; Galli, G. Strongly
Anisotropic Dielectric Relaxation of Water at the Nanoscale. J. Phys. Chem. Lett.
2013, 4, 2477–2481.

(38) Gekle, S.; Arnold, A. Comment on“Anomalous Dielectric Behavior of
Nanoconfined Electrolytic Solutions”.
Phys. Rev. Lett. 2013, 111, 089801.

(30) Schlaich, A.; Knapp, E. W.; Netz, R. R.
Water Dielectric Effects in Planar Confinement. Phys. Rev. Lett. 2016, 117,048001.

(39) Loche, P.; Wolde-Kidan, A.; Schlaich, A.;
Bonthuis, D. J.; Netz, R. R. Comment on“Hydrophobic Surface Enhances Electrostatic Interaction in Water”. Phys. Rev.
Lett. 2019, 123, 049601.

(31) Fumagalli, L.; Esfandiar, A.; Fabregas, R.; Hu, S.; Ares, P.; Janardanan, A.;
Yang, Q.; Radha, B.; Taniguchi, T.;(40) Muñoz-Santiburcio, D.; Marx, D. Chemistry in nanoconfined water. Chem. Sci.
2017, 8, 3444–3452.
12

(41) Oostenbrink, C.; Villa, A.; Mark, A. E.;
Gunsteren, W. F. V. A biomolecular force
field based on the free enthalpy of hydration and solvation: The GROMOS forcefield parameter sets 53A5 and 53A6. Journal of Computational Chemistry 2004, 25,
1656–1676.

(50) Oster, G.; Kirkwood, J. G. The Influenceof Hindered Molecular Rotation on theDielectric Constants of Water, Alcohols,
and Other Polar Liquids. J. Chem. Phys.
1943, 11, 175–178.
(51) Dresselhaus, M.;
Dresselhaus, G.;
Saito, R. Physics of carbon nanotubes.
Carbon 1995, 33, 883–891.

(42) Werder, T.; Walther, J. H.; Jaffe, R. L.;
Halicioglu, T.; Koumoutsakos, P. On
the Water-Carbon Interaction for Use
in Molecular Dynamics Simulations ofGraphite and Carbon Nanotubes. J. Phys.
Chem. B 2003, 107, 1345–1352.

(52) Loche, P.; Ayaz, C.; Schlaich, A.; Bonthuis, D. J.; Netz, R. R. Breakdown of Linear Dielectric Theory for the Interaction
between Hydrated Ions and Graphene. J.
Phys. Chem. Lett. 2018, 9, 6463–6468.

(43) Bussi, G.; Donadio, D.; Parrinello, M.
Canonical sampling through velocity
rescaling. The Journal of Chemical
Physics 2007, 126, 014101.

(53) Maggs, A. C.; Everaers, R. Simulating
Nanoscale Dielectric Response. Phys. Rev.
Lett. 2006, 96, 230603.
(54) Simoncelli, M.; Ganfoud, N.; Sene, A.;
Haefele, M.; Daffos, B.; Taberna, P.L.; Salanne, M.; Simon, P.; Rotenberg, B. Blue Energy and Desalination
with Nanoporous Carbon Electrodes: Capacitance from Molecular Simulations toContinuous Models. Phys. Rev. X 2018,8, 021024.

(44) Abraham, M. J.; Murtola, T.; Schulz, R.;
Páll, S.; Smith, J. C.; Hess, B.; Lindahl, E. GROMACS: High performancemolecular simulations through multi-level
parallelism from laptops to supercomputers. SoftwareX 2015, 1-2, 19–25.
(45) Michaud-Agrawal, N.; Denning, E. J.;
Woolf, T. B.; Beckstein, O. MDAnalysis:
A toolkit for the analysis of molecular dynamics simulations. Journal of Computational Chemistry 2011, 32, 2319–2327.
(46) Hasted, J. B. Aqueous Dielectrics (Studies
in Chemical Physics); Chapman and Hall,1974.
(47) Köfinger, J.; Dellago, C. Orientational
Dynamics and Dielectric Response of
Nanopore Water. Phys. Rev. Lett. 2009,
103, 080601.
(48) Israelachvili, J. N. Intermolecular andSurface Forces, 3rd ed.; Elsevier Science,
2010.
(49) Sedlmeier, F.; Netz, R. R. Solvation thermodynamics and heat capacity of polarand charged solutes in water. The Journal of Chemical Physics 2013, 138 .

13

Graphical TOC Entry

14