Table Of ContentOne-dimensional counterion gas between charged surfaces:
Exact results compared with weak- and strong-coupling analysis
David S. Dean,1,2 Ron R. Horgan,3,2 Ali Naji,4,2 and Rudolf Podgornik5,2
1 Universit´e de Toulouse; UPS; Laboratoire de Physique Th´eorique (IRSAMC); F-31062 Toulouse, France
2Kavli Institute of Theoretical Physics, University of California, Santa Barbara, CA 93106, USA
3DAMTP, CMS, University of Cambridge, Cambridge, CB3 0WA, UK
4Dept. of Physics, Dept. of Chemistry and Biochemistry,
and Materials Research Laboratory, University of California, Santa Barbara, CA 93106, USA
9 5Dept. of Physics, Faculty of Mathematics and Physics, University of Ljubljana,
0
and Dept. of Theoretical Physics, J. Stefan Institute, SI-1000 Ljubljana, Slovenia,
0
and Lab. of Physical and Structural Biology, National Institutes of Health, MD 20892, USA
2
(Dated: 8 January 2009)
n
a Weevaluateexactlythestatisticalintegralforaninhomogeneousone-dimensionalcounterion-only
J Coulomb gas between two charged boundaries and from this compute the effective interaction, or
disjoining pressure, betweentheboundingsurfaces. Ourexact results arecompared with thelimit-
8
ingcasesofweakandstrongcouplingwhicharethesamefor1Dand3Dsystems. Forsystemswith
] a large number of counterions it is found that the weak coupling (mean-field) approximation for
h thedisjoining pressureworksperfectly and thatfluctuationsaroundthemean-fieldin1D aremuch
c smaller than in 3D. In the case of few counterions it works less well and strong coupling approx-
e
imation performs much better as it takes into account properly the discreteness of the counterion
m
charges.
-
at PACSnumbers: 05.20.Jj,61.20.Qg
t Keywords:
s
.
t
a
m I. INTRODUCTION
-
d Electrostatic effects are predominant in a variety of soft condensed matter systems such as polymers in solutions,
n colloidal suspensions and much of the physics of membranes and films [1, 2]. The statistical mechanics of these
o
problemscannotingeneralbe solvedexplicitly;howeveranumberofapproacheshavebeendevised,startingwiththe
c
mean-field Poisson-Boltzmannequation, that allow for an approximate solution to the problem [3]. One can improve
[
the accuracy of mean-field calculations by looking at the fluctuations about the Poisson-Boltzmann solution, and in
1 a number of cases these fluctuation effects are essential to capture the quantitative behavior of the systems under
v
study. The approach of this type can be shown to be valid in the so called weak-coupling regime. More recently the
7
strong-coupling expansion, in its essence similar to the virial expansion, has been developed and successfully applied
3
to a number of interesting problems [3, 4]. Despite the success of these approaches many questions remain to be
0
1 answeredaboutthe intermediateregimeandhowthe cross-overbetweenthe weak-andstrong-couplinglimits occurs.
. Herewe showthatthe one-component1DinhomogeneousCoulombfluidcanbe analyticallysolvedandthus provides
1
an ideal test-bed to study the domains of validity of these commonly used approximationschemes.
0
9 For a system of two charged surfaces interacting in an electrolyte solution, the effective coupling strength may
0 dependonthe inter-plateseparation,andsothecalculationofthe forcebetweenthe platesrequiresanunderstanding
: of the both weak and strong coupling regimes and the transition between them. In two important papers by Lenard
v
[5] and Edwards and Lenard [6] a complete solution of the thermodynamics of the two component Coulomb gas in
i
X one dimension was given. The second of these papers was based on the fact that the underlying Sine-Gordon theory
r for this system can be mapped onto quantum mechanics in one dimension and the thermodynamics can be solved in
a
terms ofthe groundstate energy of the correspondingSchr¨odingerequation. Subsequently the system was studied in
the presence of an electric field and it was shown that it always behaves as a dielectric because of a dimerization of
the individual charge carriers [7]. When the system is not in an electroneutral state, i.e. when the surface charges
corresponding to the applied field are not exactly compensated by the counterions of the system, the physics of
the system is subtly different and confinement phenomena between oppositely charged particles appear [8, 9]. This
confinement is directly related to the fact that the system is dielectric and non-conducting. More recently, the one-
dimensional Coulomb gas thermodynamics has been studied taking into account additional hard-core and dipolar
interactions [10]. The model has also been studied in the finite size setting of a soap film type model where charge
regulation mechanism model of the surface, in the presence of salt solution (ions and counterions), is included [11].
In the present paper we study a related model to those outlined above: a one-dimensional model of counterions
confinedbetweentwoplates,carryingfixedbutnotnecessarilyequalcharges. Thissituationisquite differenttomost
of those studied above as most of them have concentrated on symmetric electrolyte systems. We place particular
2
emphasis on the computation of the pressure of the system which is the effective interaction between the two plates
modified by the presence of counterions. The effect of asymmetric charge distribution on the plates is analyzed
in detail. This system, while being somewhat idealized, is a simple model of an array of charged smectic bilayers
sandwiched between two fixed parallel charged plates. If the charge on the bilayers is uniform then the use of the
one-dimensional Coulomb interaction is justified.
Another reason for studying this simple system is to check the validity of the various approximation schemes
alluded to above. The mean-field and strong-coupling approximations in the planar geometry are independent of
the dimensionality of the problem. This is due simply to the fact the both Poisson-Boltzmann theory as well as
the strong-coupling theory are one-dimensional effective theories since the corresponding fields depend only on the
coordinate parallel to the bounding surfaces normals. However, fluctuations about the mean field, or the Poisson-
Boltmzann configuration, do depend on the dimensionality of the problem as they correspond to thermal Casimir
or zero-frequency van-der-Waals forces. This makes these model calculations particularly appealing since they can
clarify the role the fluctuations play in Coulomb systems.
The paper is organized as follows: first we use the Edwards-Lenard path integral formulation for the problem
adapted to apply to a finite electroneutral system with surface charges. The exact expression obtained for the force
between the two plates is evaluated numerically and comparedwith the results of Monte-Carlo simulation for a wide
range of cases: symmetric/asymmetric surface charges, number of counterions, etc. We find excellent agreement
between theory and simulation in all cases as we should expect since the theoretical solution is exact. In the next
section we compare the results with the Poisson-Boltzmann mean-field approximation and investigate its domain of
validity. We present two methods for calculating the effects of field fluctuations about the mean-field solution and
compare with our exact results. The strong-coupling approximation is then applied to the system and a similar
comparison with the exact results is made.
II. THE COUNTERION MODEL
We consider a system of charged particles interacting via a Coulomb potential in one dimension. It can also be
consideredas a systemof uniformly chargedinfinite two dimensionalsheets inthree dimensions,where the sheets are
perpendicular to the direction x, are uniformly charged and can only move in this direction. The Hamiltonian for
this system is given by
βe2
β = z z x x , (1)
α β α β
H − 4 | − |
α,β
X
where, e is the unit electric charge, β = 1/k T, z is the valence of particle (or sheet) α and x is its position. We
B α α
will consider systems in a canonical formulation which are overallcharge neutral, so that
z =0. (2)
α
α
X
It is convenient to rewrite the Hamiltonian as
βe2
β = z z x +x 2min(x ,x )
α β α β α β
H − 4 −
α,β
X (cid:2) (cid:3)
βe2
= z z min(x ,x ), (3)
α β α β
2
α,β
X
wherewehaveusedchargeneutralitytogofromthefirsttosecondequationabove. WemaynowwritetheHamiltonian
as an expectation value over a “Brownian motion”
2
1
β = z ψ(x ) , (4)
α α
H 2(cid:28)" α # (cid:29)
X
where ψ is Brownian motion with correlationfunction
ψ(x)ψ(y) =βe2min(x,y), (5)
h i
3
and started at the value zero at x=0, so
ψ(0)=0. (6)
The Boltzmann weight for any configuration can thus be written as
exp( β )= exp[i z ψ(x )] . (7)
α α
− H h i
α
X
In path-integral notation the measure on the Brownian motion ψ as defined by Eq. (5) is given by
∞ ψ(L) 1 L dψ(x) 2
[ψ] = dψ(L) d[ψ]exp dx [ψ],
hO i Z−∞ Zψ(0) −2e2β Z0 (cid:18) dx (cid:19) !O
(8)
where L is the overalllength of the system.
The sum over α has three distinct contributions: that from the counterions with 1 α N, for which from now
≤ ≤
on we assume have valence z = 1, and those from the surface charge at x = 0, with valence z = (N +P)/2,
α 0 0
−
and from the surface charge at x = L, with valence z = (N P)/2. The system is globally electroneutral,
N+1 N+1
− −
consistingoftwosurfacechargesandtheircounterions. Forthemoment,wewillassumethatN P areevenintegers,
±
which will be the case if the counterions originate by being released into solution from initially neutral surfaces; i.e.,
eachsurfacechargeis aninteger multiple ofthe unit ofchargeofthe counterions. The generalizationofthe approach
due to relaxing the latter assumption is given in the next section. We then have
N
N +P N P
z ψ(x )= ψ(x ) ψ(0) − ψ(L), (9)
α α i
− 2 − 2
α i=1
X X
and the Boltzmann weight for a given configuration is thus
N
N +P N P
BW( x )= exp i ψ(x ) iψ(0) iψ(L) − . (10)
i i
{ } (cid:28) i=1 − 2 − 2 !(cid:29)
X
The partition function is then obtained as
1 L N
Z = dx BW( x )
N i i
N! { }
(cid:28) Z0 i=1 (cid:29)
Y
N
1 L N +P N P
= dxexp(iψ(x)) exp iψ(0) iψ(L) − . (11)
N!(cid:28) Z0 ! (cid:18)− 2 − 2 (cid:19)(cid:29)
We now introduce an arbitrary fugacity κ having the dimensions of inverse length and write
N +P N P
κNZ = exp iψ(0) iψ(L) −
N
− 2 − 2
(cid:28) (cid:18) (cid:19)
2π dλ ∞ κN′ L N′
′
exp( iλ(N N )) dxexp(iψ(x))
×Z0 2π N′=0 − − N′! Z0 ! (cid:29)
X
2π dλ N +P N P L
= exp iψ(0) iψ(L) − iλN exp κ dxexp(iλ+iψ(x)) , (12)
(cid:28)Z0 2π (cid:18)− 2 − 2 − (cid:19) Z0 !(cid:29)
although, clearly, the value of κ should not affect the final physical results since charge neutrality will enforce the
′
constraint N =N. This is ensured by the integral over dλ. If we define a shifted Brownian motion,
φ(x)=λ+ψ(x), (13)
the starting position of the field φ at x = 0 is simply φ(0) = λ but the integration over the end point remains free
and we have
2π dφ(0) N +P N P L
Z = exp iφ(0) iφ(L) − exp κ dxexp(iφ(x)) . (14)
N
(cid:28)Z0 2π (cid:18)− 2 − 2 (cid:19) Z0 !(cid:29)
4
Theelectroneutralityconstrainthasthusbeenabsorbedintotheintegrationoverthestartingpointoftheshiftedfield
φ.
It is convenient, for what comes later, to generalize the result to the case where the counterions have valence
z =q, 1 α N. This is easily accomplished by replacing e by qe in the measure for the Brownian motion in Eq.
α
≤ ≤
(8). This follows from the forgoing analysis and rescaling the field φ(x) φ(x)/q. In path integral notation we then
→
have
2π dφ(0) ∞ N +P N P
Z = dφ(L)exp iφ(0) iφ(L) −
N
Z0 2π Z−∞ (cid:18)− 2 − 2 (cid:19)
φ(L) L 1 dφ(x) 2
d[φ]exp dx κ exp(iφ(x)) (15)
×Zφ(0) −Z0 "2q2e2β(cid:18) dx (cid:19) − #!
as our final functional integral expression for the partition function.
III. EXACT EVALUATION OF THE PARTITION FUNCTION
We use the Feynman formula, which is the Euclidean version of the equivalence between path integrals and the
propagator for the Schr¨odinger equation. We have that
∞ φ(L) L 1 dφ(x) 2 N P
dφ(L) d[φ]exp dx κ exp(iφ(x)) exp iφ(L) −
Z−∞ Zz −Z0 "2q2e2β(cid:18) dx (cid:19) − #! (cid:18)− 2 (cid:19)
N P
= exp( LHˆ)exp iz − , (16)
− − 2
(cid:18) (cid:19)
where q is the counterion valence and where Hˆ is the complex operator/Hamiltoniangiven by
q2e2β d2
Hˆ = κ exp(iz). (17)
− 2 dz2 −
The partition function can thus be expressed as
2π dz N +P N P
κNZ (L)= exp iz exp( LHˆ) exp iz − . (18)
N
2π − 2 − − 2
Z0 (cid:18) (cid:19) (cid:18) (cid:19)
We should note that the operator Hˆ is complex in our problem but in the case of symmetric electrolytes [6] it is
real. Becauseofthis we analysethe problemvia Fourieranalysisratherthanonthe basis ofthe eigenfunctions ofthe
operator Hˆ.
We define integersM =(N+P)/2 and M =(N P)/2 andthe surface chargesare M qe and M qe atx=0
1 2 1 2
− − −
and x = L, respectively, with M +M = N. A quick preliminary check of our approach is to consider the perfect
1 2
gas limit where e = 0. In this case we find trivially that
LM1+M2
κNZ (L)=κN , (19)
N
(M +M )!
1 2
and thus
LN
Z (L)= , (20)
N
N!
which is the perfect gas result. In order to proceed further we consider the evaluation of
f(z;L)=exp( LHˆ)f(z;0). (21)
−
We can write [11]
f(z)= a(n,L)exp(inz), (22)
n
X
5
and the induced evolution equation for the Fourier coefficients is
da(n,L) n2q2e2β
= a(n,L)+κa(n 1,L). (23)
dL − 2 −
In our problem we have f(z,M ;0) = exp( izM ) and thus the initial condition on the Fourier coefficients
2 2
a(n,M2;0)=δn,−M2. We define an =κnbn and−find that
db(n,M ;L) n2q2e2β
2
= b(n,M ;L)+b(n 1,M ;L), (24)
2 2
dL − 2 −
withtheinitialconditionb(n,M2;L)=κnδn,−M2 andfromEq. (18)wefindthatwecanexpressthepartitionfunction
for system of unit charges M and M on the left and right boundaries respectively as
1 2
Z (L)=b(M ,M ;L), (25)
M1,M2 1 2
but with the initial conditionbn(0,M2)=δn,−M2; as we expect, the fugacity κ, which wasintroducedon dimensional
grounds, does not enter in the final physical result. The evolution Eq. (24) can be written as
p = 1 db(M1,M2;L) = M12q2e2 + ZM1−1,M2 , (26)
M1,M2 βb(M ,M ;L) dL − 2 βZ
1 2 M1,M2
where p is the pressure of the system. Clearly, the last term makes a positive contribution to the pressure as
M1,M2
it is the ratio of two partition functions for systems with different surface charges. Equation (26) has a nice physical
interpretation. The secondterm is,up to the factorof β, the averagedensity of counterionsatz =0 because it is the
restrictedpartition function where at leastone particle is on the surface z =0 (thus giving a reductionin the surface
charge to M 1) normalized by the full partition function. We can thus write
−
σ2 1
p = 1 + ρ(0), (27)
M1,M2 − 2 β
where ρ(x) = N δ(x x ) is the average value of the counterion density at the surface z = 0 and σ = M qe
h i=1 − i i 1 − 1
the corresponding surface charge. However, this is simply the contact-value theorem for electrostatic systems known
to be exactinaPnydimension[10,12,13, 14]. Infactthe contact-valuetheoremcanbe demonstratedvia anextension
of the path integral methods used here to higher dimensional systems with (parallel) planar geometries [15]. The
average counterion density at the point x can be shown [6, 11] to be equal to
ρ(x) = exp(iφ(x))
h i
1 2π dz
= exp( iM z)exp( xHˆ)exp(iz)exp( (L x)Hˆ)exp( iM z).
1 2
Z 2π − − − − −
M1,M2 Z0
If we set x=0 in this formula we obtain
ρ(0)= ZM1−1,M2 , (28)
Z
M1,M2
in agreement with the previous discussion.
A. The general formalism
The charges on the surfaces at z = 0 and L can be denoted σ and σ respectively. Without loss of generality we
1 2
may assume that q >0 and σ σ , which allows σ and σ to be of opposite sign, and we can define the charge
2 1 1 2
| |≥ | |
asymmetryparameterζ =σ /σ with 1 ζ 1. Allothercasescanbe mappedontothis intervalwithappropriate
2 1
− ≤ ≤
rescaling of the parameters. Note that the values ζ =1 and ζ = 1 representspecial cases of symmetric system with
−
σ = σ , and antisymmetric system with σ = σ and no counterions between surfaces (which thus reduces to the
1 2 1 2
−
trivial case of a planar capacitor).
For the above analysis we see that
N P
ζ = − , (29)
N +P
6
σ = M qe and σ = M qe. The independent parameters are conveniently chosen as σ ,ζ,N. With M and M
1 1 2 2 1 1 2
− −
integers, the analysis so far allows only for particular discrete values of ζ. However, the approach can be extended
to accommodate all values of ζ and we state the generalization of the algorithm here. Charge neutrality determines
that the valence q of the counterion satisfies
σ +σ σ
1 2 1
q = = (1+ζ). (30)
− Ne −Ne
We are allowing that q need not be an integer.
In additionto the Fourier coefficients b(n,M ;L)we introduce Fouriercoefficients c(n,M ;L). Equation(24) gives
2 1
the evolution of the b(n,M ;L) from z =0, the surface with assigned charge σ . The coefficients c(n,M ;L) obey a
2 1 1
similarevolutionequationbutnowwiththe surfacechargeassignmentreversed;the chargeatz =0is nowσ . These
2
twosetsofcoefficientscorrespondtocomplementaryapproaches;theset b describesevolutionfromthesurfacewith
{ }
charge σ to that with chargeσ , and vice-versa for the set c .
1 2
{ }
We define
1
α= , M =Int(αN), η =αN M , η =1 η , M =N M 1. (31)
1 1 1 2 1 2 1
1+ζ − − − −
The evolution equations now become
db(n,M ;L) (n η )2
2 = − 2 βq2e2b(n,M ;L)+b(n 1,M ;L), (32)
2 2
dL 2 −
dc(n,M ;L) (n η )2
1 = − 1 βq2e2c(n,M ;L)+c(n 1,M ;L). (33)
1 1
dL 2 −
The initial conditions are b(n,M2;L=0)=δn,−M2, c(n,M1;L=0)=δn,−M1.
Then the partition function is given by alternative formulas Z(σ ,ζ,N;L)=b(M +1,M ;L)=c(M +1,M ;L),
1 1 2 2 1
and the pressure is given equivalently by
(M +1 η )2 b(M ,M ;L)
1 − 2 q2e2 + 1 2
− 2 Z(σ ,ζ,N;L)
1
p(L)= (34)
(M2+1−η1)2q2e2 + c(M2,M1;L) .
− 2 Z(σ ,ζ,N;L)
1
The partition function Z(σ1,ζ,N;L)is also given by
M2+1
Z(σ ,ζ,N;L)= b( j+1,M ;L x)c(j,M ;x), (35)
1 2 1
− −
j=X−M1
for any x, 0 x L. This acts as a check on the numerics.
≤ ≤
The counterion number density ρ(x;σ ,ζ,N,L) is given by
1
1 M2+1
ρ(x;σ ,ζ,N,L)= b( j,M ;L x)c(j,M ;x), (36)
1 2 1
Z(σ ,ζ,N;L) − −
1 j=X−M1
and another check on the numerics is that
dxρ(x;σ ,ζ,N,L) = N, (37)
1
Z
i.e. the system is electroneutral.
IV. APPROXIMATE EVALUATIONS OF THE PARTITION FUNCTION
A traditional approach to the study of charged (bio)colloidal systems is the mean-field Poisson-Boltzmann (PB)
formalismwhich is applicable for weak surface charges,lowcounter-ionvalency and high temperature [16]. The limi-
tations of this approach are evident when applied to highly-charged systems where counterion-mediated interactions
7
between charged bodies start to deviate substantially from the accepted mean-field wisdom [3, 4]. One of the recent
fundamental advances in this field has been the systematization of these non-PB effects based on the notions of weak
and strong coupling approximations. The latter approachhas been pioneered by Rouzina and Bloomfield [17], elabo-
rated later by Shklovskii et al. [18], Levin et al. [19], and brought into final form by Netz et al. [3, 4, 20, 21]. These
two approximations allow for an explicit and exact treatment of charged systems at two disjoint limiting conditions
whereas the parameter space in-between can be analyzed only approximately [20, 21, 22, 23, 24, 25] and is mostly
accessible solely via computer simulations [3, 4, 20, 21, 24, 25, 26, 27, 28, 29, 30, 31].
Both the weak- and the strong-coupling approximations are based on a functional integral or field-theoretic rep-
resentation [32, 33] of the grand canonical partition function for a system composed of fixed surface charges with
intervening mobile counterions, and depend on the value of a single dimensionless coupling parameter Ξ [20, 21]. In
three dimensions, Ξ is proportional to the ratio between two relevant length scales, namely, the Bjerrum length and
the Gouy-Chapman length. The Bjerrum length, defined as ℓ = e2/(4πεε k T), is the distance at which two unit
B 0 B
chargesinteractwiththermalenergyk T (inwateratroomtemperature,one hasℓ 0.7nm). Ifthe chargevalency
B B
of the counterionsis q then the correspondinglength scales as q2ℓ . Similarly, the Go≃uy-Chapmanlength, defined as
B
µ = e/(2πqℓ σ ), is the distance at which a counterion interacts with a macromolecular surface (of surface charge
B
| |
density σ) with an energy equal to k T. The 3D electrostatic coupling parameter measures the competition between
B
ion-ion and ion-surface interactions and is given by [20, 21]
Ξ=Ξ q2ℓ /µ=2πq3ℓ2 σ /e. (38)
3D ≡ B B| |
Physically, the weak-coupling (WC) regime Ξ 1 (appropriate for low valency counterions and/or weakly charged
≪
surfaces), is characterizedby the fact that the width, µ, of the counterion layer near the surfaces is much larger than
the separationbetween two neighboring counterions in solution, and thus the counterion layer behaves basically as a
three-dimensionalgas. Eachcounterioninthiscaseinteractswithmanyothersandthecollectivemean-fieldapproach
of the Poisson-Boltzmann(PB) type is completely justified.
On the other hand in the strong-coupling (SC) regime Ξ 1 (appropriate for high valency counterions and/or
≫
highly charged surfaces), the mean distance between counterions, a⊥ = σ /qe, is much larger than the layer width
| |
(i.e., a⊥/µ √Ξ 1), indicating that the counterions are highly localized laterally and form a strongly correlated
∼ ≫
quasi-two-dimensional layer next to a charged surface. In this case, the weak-coupling approach breaks down due to
strong counterion-surface and counterion-counterion correlations. Since counterions can move almost independently
from the others along the direction perpendicular to the surface, the collective many-body effects that enable a
mean-field description are absent, necessitating a complementary SC description [20, 21].
Formally, the WC limit can be identified with the saddle-point approximation of the field theoretic representation
ofthegrandcanonicalpartitionfunction,andreducestothemean-fieldPBtheoryatthe lowestorderforΞ 0. The
→
quadratic fluctuations around the mean field provide a second-order correction to the mean-field solution for small
Ξ < 1 [32, 33, 34, 35, 36, 37, 38, 39]. The SC approximation has no PB-like collective mean-field [20, 21] since it is
formally equivalent to a single particle description obtained from a systematic 1/Ξ expansion in the limit Ξ ,
→ ∞
and corresponds to two lowest order terms in the virial expansion of the grand canonical partition function. The
consequences and the formalism of these two limits of the Coulomb fluid description have been explored widely and
in detail (for reviews, see Refs. [3, 4]).
The concept of weak- and strong-coupling electrostatics can be easily generalized to other dimensions. In one
dimension (1D), the Coulomb interaction between two unit charges may be written as v (x) = xk T/ℓ , where
1D B B
ℓ =2k T/e2mayberegardedasa1DBjerrumlength. Likewise,theinteractionstrengthwithth−e|bo|undarycharges
B B
inanasymmetricsystemmaybecharacterizedbythe1DGouy-Chapmanlengthsµ =ℓ e/q σ andµ =ℓ e/q σ .
1 B 1 2 B 2
| | | |
Wemayproceedbyusingonlyµ µ =ℓ e/q σ inwhatfollowssinceforanygivenasymmetryparameterζ =σ /σ ,
1 B 1 2 1
≡ | |
we have µ =µ/ζ . The ratio
2
| |
Ξ q2µ/ℓ (39)
1D B
≡
is then the corresponding electrostatic coupling parameter in 1D. (On a more formal level, this definition for the
coupling parameter may be established by looking at the field-theoretical representation(15) and rescaling the coor-
dinates as x x /µ .) Note that in 1D the coupling parameter is related to the asymmetry parameter by virtue of
i i 1
→
the electroneutrality condition
Nqe= (σ +σ ), (40)
1 2
−
where N is the number of counterions of valency q, as
ζ+1
Ξ = . (41)
1D
N
8
Thus, it may be expected that in the limit N , contrary to the 3D case, the mean-field theory corresponding
→∞
to this system becomes exact! On the other hand, the SC description may be expected to follow simply for N = 1.
These limiting caseswill be discussedfurther inthe forthcomingsectionsandwill be comparedwith the exactresults
and Monte-Carlo simulations. One should bear in mind that the 1D system considered here corresponds to a 3D
system of mobile charged plates (membranes) confined between two fixed planar charged walls and thus, such a 1D
system with even a single “counterion” would have a completely meaningful thermodynamic behavior.
We note also that both the PB theory (Ξ 0) and the SC theory (Ξ ) for uniformly charged plates are one-
→ →∞
dimensionaltheoriesandshouldremainvalidin3Daswellasinthe1Dcasethatwearestudyinghere. However,even
though the mean-field solution depends only on the transversecoordinate,the field fluctuations about the mean-field
solution will depend on all coordinates and so it is clear that fluctuations corrections to the mean-field contribution
depend on the dimensionality of the system: they will be different for a 1D system than for a 3D system.
A. Weak-coupling limit: Poisson-Boltzmann theory
The expression for the partition function Z in Eq. (15) is up to multiplicative constants given by the functional
N
integral
Z = d[φ]exp( S[φ]) (42)
N
−
Z
where the action S is given by
L 1 dφ 2 σ σ
1 2
S[φ]= dx κexp(iφ) iφ δ(z) iφ δ(z L) (43)
Z0 "2q2e2β (cid:18)dx(cid:19) − − qe − qe − #
In the weak-coupling regime, the leading contribution to the partition function comes from the saddle-point config-
uration, φ (x), of the action in Eq. (43) [32, 33], where the field φ (x) = iψ (x) turns out to be imaginary and is
0 0 0
proportionalto the mean-field electrostatic potential. The saddle-point configurationcan be straightforwardlytrans-
latedintoasolutionofthe PBequationandcorrespondstoanexactasymptoticresultinthe limitΞ 0[20,21,34].
→
The PB equation for the potential, ψ (z), can be written as [2, 16]
0
d2ψ (x)
0 = q2e2βκe−ψ0(x), (44)
dx2 −
with boundary conditions
dψ 2
0
= σ βqe= ,
1
dx − µ
(cid:12)0
(cid:12)
dψ0(cid:12) 2ζ
(cid:12) =σ2βqe= (45)
dx − µ
(cid:12)L
(cid:12)
(cid:12)
stemming from the electroneutrality of the system. (cid:12)
Integration of the PB equation gives rise to the first integral of the system of the form
2
1 dψ
0
βp = +ρ (x), (46)
0 −2q2e2β dx 0
(cid:18) (cid:19)
where the constant p is nothing but the mean-field PB pressure acting between the bounding surfaces [2] and
0
ρ (x)=κe−ψ0(x), (47)
0
isthePBnumberdensityprofileofcounterionsbetweenthesurfaces. AsremarkedearlierinsectionIII,thechoicefor
κ is arbitrary since it corresponds to a choice of origin for the potential ψ(x). In section III we couched the Fourier
solution in terms of Fourier coefficients a =κnb with the choice κ=1, and so we adopt this choice here, too.
n n
The natureofthe solutionψ (x)obviouslycruciallydepends onthe signofthe pressurep [40,41]. Differentforms
0 0
are obtained for positive and negative pressures, corresponding to repulsion and attraction between the bounding
surfaces respectively and were derived previously [42].
9
Guided by the contact-value theorem Eq. (27), which we emphasize is an exact result, we express the Poisson-
Boltzmann, or mean-field, pressure in rescaled units. The repulsive PB pressure is given by
2p βp ℓ
p˜ = 0 = 0 B =α˜2. (48)
0 σ2 (σ /e)2
1 1
and all the other quantities have been defined above. Furthermore, α˜ =αµ and its value is given by the solution of
α(ζ +1)µ
tan(αL)= . (49)
α2µ2 ζ
−
For ζ <0, and only then, can the PB pressure be attractive. The attractive PB pressure is given by
2p βp ℓ
p˜ = 0 = 0 B = α˜2. (50)
0 σ2 (σ /e)2 −
1 1
where α˜ is now given as a solution of
ζ+µ2α2
coth(αL)= . (51)
−µα(1+ζ)
Thelimitingcaseofzeropressurecanbe obtainedstraightforwardlyfromeitherlimit. Theinteractionpressurecan
be obviously computed for any value of the asymmetry parameter, but note that any value of ζ can be mapped onto
the interval 1<ζ <1 and the pressure only need to be evaluated in that interval of ζ values. Note again that the
−
PB interaction pressure is the same for a 1D as well as for a 3D system.
B. Weak-coupling limit: fluctuations
Fluctuations around the mean field depend on the dimensionality of the system and are different for a 1D than for
a 3D system. In order to proceed, one needs to evaluate the appropriate Hessian of the field action in the partition
function and study its fluctuation spectrum (see Refs. [32, 33, 43] for more details). The Hessian of the field about
the mean-field value is given by
δ2S 1 d2
=H(x,x′)= +βe2q2ρ (x) δ(x x′), (52)
δφ(x)φ(y)|φ=iψ0 βe2q2 −dx2 0 −
(cid:20) (cid:21)
ρ (x)=exp( ψ (x)) is the mean-field PB counterion density. Hence,
0 0
−
2α2
p >0,
cos2α(x x ) 0
βe2q2ρ (x)= − 0 (53)
0 2α2 p <0.
sinh2α(x x ) 0
0
The corresponding correction, , to the free energy of the syst−em comes from the functional integral
2
F
ψ(L)=y
(L)= dxdy d[ψ]exp( s[ψ]), (54)
A −
Z Zψ(0)=x
where
1
′ ′ ′
s[ψ]= ψ(x)H(x,x)ψ(x)dxdx, (55)
2
ZZ
so that the corresponding fluctuation part of the free energy can be obtained in the form
β = ln (L). (56)
2
F − A
The functional integral Eq. (54) can be evaluated exactly, in two different ways. The first method is based on
the use of the argument principle [36] converting the discrete sum of eigenvalues of the Hessian operator into the
10
logarithm of the secular determinant of the same operator. The second method is based on the Pauli-van Vleck
λ
D
approach to calculating the functional integral of a general harmonic kernel. Both results are the same.
In the first approachthe trace-logof the Hessiancan be written equivalently in the form that was derivedfor a 3D
case [32, 33, 36] but can be used in a trivially modified form also for the 1D case under consideration. It gives
β = 1Tr lnH(x,x′)= 1lnD1, (57)
F2 2 2
0
D
where istheseculardeterminant,i.e. thedeterminantofthecoefficientscorrespondingtotheappropriateboundary
λ
D
condition in the solution of the eigenvalue equation, that can be derived from solutions of
d2
λβe2q2ρ (x) f (x)=0. (58)
dx2 − 0 λ
(cid:16) (cid:17)
One can now write (L) for the quotient ( / )−1/2, since the secular determinant depends explicitly on the value
1 0
A D D
of the inter-surface spacing, L. Using the fact [32, 33, 43] that the two linearly independent solutions of Eq. (58) for
λ=1 in the repulsive regime are
f(1)(x)=tanαx and f(2)(x)=1+αxtanαx (59)
the secular determinant for the symmetric case ζ =1 comes out as
L L L L
(L)=αsec2 α tan α +α sec2 α , (60)
1
D 2 2 2 2
(cid:18) (cid:19)(cid:20) (cid:18) (cid:19) (cid:18) (cid:19)(cid:21)
and the corresponding free energy contribution from the quadratic fluctuations is thus given by
1 L L L L
β = ln αsec2 α tan α +α sec2 α . (61)
2
F 2 2 2 2 2
(cid:26) (cid:18) (cid:19)(cid:20) (cid:18) (cid:19) (cid:18) (cid:19)(cid:21)(cid:27)
The fluctuation free energy can be regularized so that all irrelevant constants, i.e. all the terms not depending on
the separation between the bounding surfaces, are dropped, amounting to a rescaling (L) (L) (L ).
2 2 2
F →F −F →∞
This corresponds to a subtraction of the part of the free energy for two separate boundaries at infinite separation
from the total free energy.
In the second method we compute the fluctuations about the mean-field path using the Pauli-van Vleck approach
[44, 45, 46]. As before the thickness of the film is L we take the leftmost and rightmost points of the film to be at
x=0andx=Lrespectively. Clearlythe actionofthe classicalpathminimizings inEq. (54)isaquadraticfunction
of the initial and final points of the fluctuating field ψ. The generalized Pauli-van Vleck formula tells us that
1
ψ(L)=y 1 ∂s [x,y] 2
c
d[ψ]exp( s[ψ])= exp( s [x,y]), (62)
c
− −2π ∂x∂y −
Zψ(0)=x (cid:18) (cid:19)
where s is the classical action minimizing s. As the action is quadratic we may write
c
1
s [x,y]= f(L)x2+g(L)y2 2h(L)xy . (63)
c
2 −
(cid:2) (cid:3)
In the case of a symmetric charge distribution we have f =g and thus the fluctuation term is given by
1
2πh(L) 2
(L)= . (64)
A f2(L) h2(L)
(cid:18) − (cid:19)
Solving the equations of motion we find that
αsec2(αL/2) α[tan(αL/2)+(αL/2)sec2(αL/2)]
βq2e2f(L) = βq2e2g(L)= + ,
2tan(αL/2) 2[1+(αL/2)tan(αL/2)]
α
βq2e2h(L) = . (65)
2tan(αL/2)[1+(αL/2)tan(αL/2)]
Putting all this together then yields
1
π 2
(L)= (66)
A αsec2(αL/2)[tan(αL/2)+(αL/2)sec2(αL/2)]
(cid:18) (cid:19)
