Table Of ContentElectrostatics on the sphere with applications to Monte Carlo
simulations of two dimensional polar fluids.
Jean-Michel Caillol∗
University of Paris-Sud, CNRS, LPT,
5
1 UMR 8627, Orsay, F-91405, France
0
2 (Dated: January 23, 2015)
n
a Abstract
J
2 We presenttwo methodsfor solving theelectrostatics of pointcharges andmultipoles on thesur-
2
] face of a sphere, i.e. in thespace S2, with applications tonumerical simulations of two-dimensional
t
f polar fluids.
o
s
. In the first approach, point charges are associated with uniform neutralizing backgrounds to
t
a
m form neutral pseudo-charges, while, in the second, one instead considers bi-charges, i.e. dumbells
-
d ofantipodalpointchargesofoppositesigns. Weestablishtheexpressionsoftheelectricpotentialsof
n
o
pseudo- and bi-charges as isotropic solutions of the Laplace-Beltrami equation in . A multipolar
c S2
[
expansion of pseudo- and bi-charge potentials leads to the electric potentials of mono- and bi-
1
v multipolesrespectively. Thesepotentials constitutenon-isotropicsolutionsoftheLaplace-Beltrami
8
3 equationthegeneralsolutionofwhichinsphericalcoordinatesisrecastunderanewappealingform.
5
5
We then focus on the case of mono- and bi-dipoles and build the theory of dielectric media in
0
.
1 . We notably obtain the expression of the static dielectric constant of a uniform isotropic polar
2
0 S
5 fluid living in in term of the polarization fluctuations of subdomains of . We also derive
2 2
1 S S
:
v the long range behavior of the equilibrium pair correlation function under the assumption that it
i
X
is governed by macroscopic electrostatics. These theoretical developments find their application
r
a
in Monte Carlo simulations of the 2D fluid of dipolar hard spheres. Some preliminary numerical
experiments are discussed with a special emphasis on finite size effects, a careful study of the
thermodynamic limit, and a check of the theoretical predictions for the asymptotic behavior of the
pair correlation function.
Keywords: Laplace-Beltramiequationonthesphere;Two-dimensionalpolarfluids;MonteCarlosimulations.
∗
Electronic address: [email protected]
1
I. INTRODUCTION
The idea of using the two dimensional (2D) surface of a sphere, i.e. the space , to
2
S
perform numerical simulations of a 2D fluid phase can be tracked back to a paper by J.-P.
Hansen et al. devoted to a study of the electron gas at the surface of liquid Helium [1].
Subsequently, the same idea was used to study the crystallization of the 2D one-component
plasma (with logr interactions) [2], to establish the phase diagram of the 2D Coulomb
∝
gas [3], and to determine some thermodynamic and structural properties of the 2D polar
fluid (with 1/r2 interactions) in its liquid phase [4, 5]. The generalizations to 3D systems,
∝
implying the use of the surface of a 4D hypersphere, i.e. the space , is due to Caillol and
3
S
Levesque [6]. Some improvements on this early work, mostly applications concerning 3D
polar fluids, were published recently [7].
Several simple ideas can be put forward to justify the use of a hyperspherical geometry
in numerical simulations of plasmas and Coulombic fluids (i.e. fluids made of ions and/or
dipoles) :
(i) The n-dimensional non-Euclidian space (in practice n = 2,3), albeit finite, is
n
• S
homogeneous and isotropic, in the sense that it is invariant under the group (n+1)
O
of the (n + 1)D rotations of the Euclidian space E ; it is thus well suited for the
n+1
simulation of a fluid phase (in the bulk).
(ii) The laws of electrostatics can easily be obtained in and the Green’s function of
n
• S
Laplace-Beltrami equation (i.e. Coulomb potential) is known [8]. It has a very simple
analytical expression, tailor-made for numerical evaluations.
(iii) The inclusion of charged or uncharged walls can easily be done in order to study
•
the structure of liquids, plasmas, or collo¨ıds at interfaces [6, 9–12]
The second point (ii) was partly overlooked by the authors of Ref. [4, 5] who used an em-
pirical 2D dipole/dipole interaction, unfortunately not deduced from a solution of Laplace-
Beltrami equation in , which casts some doubts on the validity of their results concerning
2
S
the 2D Stockmayer fluid, notably those concerning its dielectric properties.
Some years after we are now in position to propose in this paper two possible, both
rigorously correct, dipole/dipole pair-potentials and to use both of them into actual numer-
ical simulations of a 2D dipolar hard sphere (DHS) fluid. Here, not only we extend to
2
S
2
the recent developments on 3D polar fluids in discussed in Ref. [7], but we also present
3
S
additional new results on multipolar expansions in , not yet available in , in Sec. II.
2 3
S S
This article is thus a contribution to the physics of 2D fluids, collo¨ıds or plasmas, the
atoms or molecules of which interact via electrostatic pair potentials derived from a solution
of the 2D Laplace equation. We emphasize that the systems that we consider cannot be seen
asthinlayers ofreal3D systems. Inthiscase, theelectrostaticinteractions shouldbederived
from the solutions of the 3D Laplace equation. Simulations of such 3D systems constrained
to live in a 2D geometry are nevertheless also possible in spherical geometries. In that case
one should consider a 3D fluid of made of particles interacting by pair potentials deduced
3
S
from solutions of Laplace-Beltrami equation in and constrained to stay on the equator of
3
S
the hypersphere. This geometrical locus indeed identifies with the sphere . This approach
2
S
was used for the 3D restricted primitive model of electrolytes in Ref. [13].
In this work our aim is to demonstate the validity of the sphere to perform Monte
2
S
Carlo (MC) simulations of a 2D polar fluid. We note that, although 2D dipolar fluids do
not exist per se in nature, the model could be used via various mappings for applications
as, recently, for the hydrodynamics of two-dimensional microfluids of droplets [14].
Our paper is organized as follows. After this introduction (Sec. I) we discuss in depth
the electrostatics of the distributions of charges in and their multipole expansions. The
2
S
material reported in Sec. II is an intricate mixture of old and new stuff. The underlying idea
is the remark of Landau and Lifchitz in Ref. [15] that, in a finite space such as , Laplace-
2
S
Beltrami equation admits solutions if and only if the total electric charge of the space is
equal to zero. Therefore, in , the building brick of electrostatics cannot be a single point
2
S
charge q as in the ordinary 2D Euclidian space E . We are led instead to consider a pseudo-
2
charge [6], a neologism denoting the association of a point charge and a uniform neutralizing
background of opposite charge. Alternatively we can consider a bi-charge, i.e. a neutral
dumbell made of two antipodal charges of opposite signs +q and q as in Ref. [16]. Both
−
approaches yield independent electrostatics theories which are sketched in Sec. II and III.
The case of dipoles is then examined more specifically in Sec. III. We need make a
distinction between mono-dipoles, which are obtained as the leading term of a multipole
expansion of a neutral set of pseudo-charges, and bi-dipoles, obtained in a similar way from
bi-charges. As a consequence two possible microscopic models of dielectric media can be
worked out, those made of mono-dipoles living on the entire sphere and those made of
2
S
3
bi-dipoles living on the northern hemisphere +.
S2
These two types of models are then studied in the framework of Fulton’s theory [17–20]
in Sec. IV. Fulton’s theory realizes a harmonious synthesis of linear response theory and the
macroscopic theory of dielectric media. It allows us to obtain
a family of formula for the dielectric constant ǫ of the fluid, well adapted for its
•
determination in MC simulations.
an expression for the asymptotic pair correlation in under the assumption that it
2
• S
should be dictated by the laws of macroscopic electrostatics.
Monte-Carlo (MC) simulations of the 2D DHS fluid in the isotropic fluid phase are then
discussed in Sec. V. We have chosen to report MC data only for two thermodynamic states,
both in the isotropic fluid phase and extensive MC simulations of the 2D DHS fluid will
be published elsewhere. These states could serve as benchmarks for future numeric studies.
One of these states was studied many years ago in Ref. [21] by means of MC simulations
within standard periodic boundary conditions (see also [22]). Comparisons of our results
with those obtained in this pioneer work are OK. A carefull finite size scaling study of our
data yields results of a high precision, notably for the energy, probably unattainable by
standard simulation methods.
Conclusions are drawn in Sec (VI).
II. ELECTROSTATICS OF 2D CHARGES AND MULTIPOLES
A. The Plane
Let us first recall that the electrostatic potential V (ρ) at a source-free point ρ = (x,y)
E2
of the Euclidian plane E satisfies the 2D Laplace equation ∆ V = 0, which in polar
2 E2 E2
coordinates (ρ,ϕ) reads
1 ∂ ∂V 1 ∂2V
∆ V = ρ + = 0 . (1)
E2 E2 ρ∂ρ ∂ρ ρ2 ∂ϕ2
(cid:18) (cid:19)
The general solution of Eq. (1) is [23, 24]
∞ ∞
V (ρ,ϕ) = a +b logρ+ a ρncos(nϕ+α )+ b ρ−ncos(nϕ+β ) , (2)
E2 0 0 n n n n
n=1 n=1
X X
4
where a , b , α , and β are arbitrary constants. The terms of the r.h.s. of (2) which are
n n n n
singular at ρ = 0 may be interpreted as the potentials created by point multipoles located
at the origin O, while those singular at ρ = as the potentials of point multipoles at
∞
infinity. In E the potential of a point charge q located at O is qlogρ up to an additional
2
−
constant; it is the Green’s function of Eq. (1) in E (without boundaries) and it satisfies the
2
2D Poisson’s equation
∆ ( logρ) = 2πδ (M,M ) 2πδ(x)δ(y) , (3)
E2 − − E2 0 ≡ −
where δ(x) denotes the 1D Dirac’s distribution.
′
Finally, the Green’s function log( ρ ρ ) can be expanded in polar coordinates as [23]
− | − |
∞ n
1 ρ
log( ρ ρ′ ) = logρ + < cos n ϕ ϕ′ , (4)
<
− | − | − n ρ −
Xn=1 (cid:18) >(cid:19) (cid:16) (cid:16) (cid:17)(cid:17)
′ ′
where ρ = inf(ρ,ρ ) and ρ = sup(ρ,ρ ).
< >
Eq. (4) serves as a starting point for the 2D multipolar expansion of Ref. [24]. Let us
consider N charges q of polar coordinates (ρ ,ϕ ) and a point M (coordinates (ρ,ϕ)) of the
i i i
plane E such that ρ > ρ , i. Making use of (4) one finds
2 i
∀
N
V (M) = q log( ρ ρ )
E2 − i | − i|
i=1
X
∞
1 1 1
= V Q logρ+ Q exp( inνϕ) , (5)
0 − 2 0ν 2 nρn nν −
ν n=1 ν
X XX
where V is some unessential constant. In (5) the index ν assumes the values 1 and +1
0
−
only. The complex multipolar moments are defined as
N
Q = q ρnexp(inνϕ ) . (6)
nν i i i
i=1
X
Since Q = Q∗ there are at most two independent multipole components at a given
n,−ν nν
order n. The potential of a 2D multipole of order n therefore decays as 1/ρn. This remark
justifies a posteriori the interpretation that we gave for the terms of the r.h.s. of Eq. (2).
B. The Sphere
Here, we examine how the basic laws of 2D Euclidian electrostatics swiftly evoked in
Sec. IIA are changed when the system is wrapped on the surface of a sphere. Let us
5
first introduce some notations and recall some elementary mathematics. We denote by
(O,R) the sphere of center O and radius R of the usual 3D geometry. It is a (Riemannian)
2
S
manifold of the 3D Euclidean space E that we identify with R3. When we deal with the
3
sphere of unit radius we adopt the uncluttered notation (O,R = 1). Let M be the
2 2
S ≡ S
running point of 2(O,R), we define the unit vector z 2 as OM = Rz. The spherical
S ∈ S
coordinates are defined as usual : z = (sinθcosϕ,sinθsinϕ,cosθ)T with 0 θ π and
≤ ≤
0 ϕ < 2π. The differential vector dz = dθe + sinθdϕe allows us to obtain two
θ ϕ
≤
orthogonal unit vectors (e ,e ) that span the plane (M) tangent to (O,R) at point M.
θ ϕ 2 2
T S
Recall that e = (cosθcosϕ,cosθsinϕ, sinθ)T and e = ( sinϕ,cosϕ,0)T. In addition
θ ϕ
− −
and quite specifically since this notion cannot be generalized to higher dimensions, one has
e = z e where the symbol denotes the 3D vectorial product. Finally the infinitesimal
ϕ θ
× ×
surface element is dS = R2dΩ where the infinitesimal solid angle dΩ = sinθdθdϕ in spherical
coordinates.
In the sequel we will make a repeated use of the unit dyadic tensor U (z) = e e +e e
S2 θ θ ϕ ϕ
of the plane (M). Note that the unit dyadic tensor of Euclidian space E is given by
2 3
T
U = U (z) + zz. It is a constant tensor independent of point M. These admittedly
E3 S2
old-fashioned objects however allow an easy definition of the gradient in (O,R), or first
2
S
differential Beltrami operator, as
= U (z) ,
∇S2(O,R) S2 ·∇E3
where is the usual Euclidian gradient operator of E and the dot in the r.h.s. denotes
∇E3 3
the 3D tensorial contraction. Note that, obviously, = /R.
∇S2(O,R) ∇S2
The Laplace-Beltrami operator (or second differential Beltrami operator) is defined in a
similar way as the restriction of the 3D Euclidian Laplacian ∆ to the surface of the sphere
E3
[25].
∂2 2 ∂
∆ = ∆ . (7)
S2(0,R) E3 − ∂R2 − R∂R
We have the scaling relation ∆ ∆ /R2 and, on the unit sphere, in spherical coor-
S2(0,R) ≡ S2
dinates
1 ∂ ∂ 1 ∂2
∆ = sinθ + . (8)
S2 sinθ∂θ ∂θ sin2θ∂ϕ2
(cid:18) (cid:19)
∆ identifies with minus the squared angular momentum operator of Quantum Mechanics;
S2
therefore it is a Hermitian operator the eigenvectors of which are the 3D spherical harmonics
6
Ym(θ,ϕ) where l is a non-negative integer and the integer m satisfies l m +l.
l − ≤ ≤
Moreover ∆ Ym = l(l +1)Ym.
S2 l − l
A general solution of Laplace-Beltrami equation in , i.e. ∆ V = 0, can easily be
S2 S2 S2
obtained in spherical coordinates by expanding V (θ,ϕ) in Fourier series (see e.g. [26])
S2
m=+∞
V (θ,ϕ) = Vˆ (θ)eimϕ , (9)
S2 m
m=−∞
X
which yields for the Fourier coefficients
d d
sinθ sinθ Vˆ (θ) = m2Vˆ (θ). (10)
m m
dθ dθ
In the above equation the change of variables x = logtan θ (0 θ π) leads to the
2 ≤ ≤
elementary differential equation
d2
Vˆ (x) = m2Vˆ (x) , (11)
dx2 m m
with the solutions
ˆ
V (x) = a+bx, for m = 0 , (12a)
0
Vˆ (x) = K e±mx, for m = 0 , (12b)
m ±m
6
where a, b, and K are arbitrary complex constants. The general real solution of Laplace-
±m
Beltrami equation on the sphere can thus be finally written as
∞ n
θ θ
V (θ,ϕ) = a +b logtan + a tan cos(nϕ+α )
S2 0 0 2 n 2 n
n=1 (cid:20) (cid:21)
X
∞ n
θ
+ b cot cos(nϕ+β ) , (13)
n n
2
n=1 (cid:20) (cid:21)
X
where a , b , α , and β are arbitrary real constants. The terms of the r.h.s. of (13)
n n n n
which are singular at θ = 0 may be interpreted as the potentials created by point multipoles
located at the north pole and those singular at θ = π as the potentials of point multipoles
N
located at the south pole . Remarkably the isotropic term log(tanθ/2) in the r.h.s. is
S ∝
singular both at θ = 0 and θ = π and should identify with the potential of a unit bi-charge.
In order to check this assertion, let us first consider a bi-charge made of a charge +q
located at point M of (0,R) and its companion q located at the antipodal point M
0 2 0
S −
(with z = z ). The potential Vbi (M) at point M is the solution of Poisson’s equation
0 − 0 q,M0
∆ Vbi (M) = 2πq δ (z ,z) δ (z ,z) , (14)
S2 q,M0 − { S2 0 − S2 0 }
7
where the Dirac’s distribution in is defined as δ (z ,z) = δ(1 z z) [25]. Eq. (14) is
S2 S2 0 − 0 ·
solved by expanding both sides on spherical harmonics; after some elementary algebra one
finds
2l+1
Vbi (M) = q ′ P (cosψ )
q,M0 l(l+1) l M0M
l
X
ψ
= q qlogtan M0M , (15)
− − 2
where P (x) is a Legendre polynomial and the prime affixed to the sum in (15) means the
l
restrictionthat l isanodd, positive integer. Inthisequationwehave introduced thegeodesic
distance ψ = arccos(z z) between the two points M and M on the unit sphere .
M0M 0 · 0 S2
The first term in the r.h.s. of Eq. (13) is thus indeed the potential created at a source-free
point M by a point bi-charge located at the north pole .
N
It is the place to introduce Dirac’s function ona sphere of radius R = 1; it will bedenoted
6
δ(M ,M) = R−2δ (z ,z). Since the Laplacian also scales as R−2 with the radius of the
0 S2 0
sphere the potential Vbi (M) is independent of the latter. Thus, in (0,R), Poisson’s
q,M0 S2
equation for a bi-charge takes the form
∆ Vbi (M) = 2πq δ(M ,M) δ(M ,M) , (16)
S2(O,R) q,M0 − 0 − 0
(cid:8) (cid:9)
The electric field given by
1
Ebi (M) = Vbi (M)
q,M0 −R∇S2 q,M0
q 1
= t (M) , (17)
Rsinψ M0M
M0M
where t (M) = cotψ z z /sinψ denotes the unit vector, tangent to the
M0M M0M − 0 M0M
geodesics M M at point M, and orientated from point M towards point M [4, 5]. Note that
0 0
this vector differs from t (M ) = cotψ z + z/sinψ which is the unit vector,
M0M 0 − M0M 0 M0M
tangent to the geodesics M M at point M (also orientated from M to point M). One
0 0 0
checks that the electric field Ebi (M) satisfies to Gauss theorem [4, 5].
q,M0
Let us consider now a pseudo-charge made of a point charge +q located at at point M
0
of and a uniform neutralizing background of charge density q/(4π). The potential
2
S −
Vps (M) at point M is the solution of Poisson’s equation
q,M0
1
∆ Vps (M) = 2πq δ (z ,z) . (18)
S2 q,M0 − { S2 0 − 4π}
8
Eq. (18) can also be solved by expanding both sides on the complete basis set of spherical
harmonics with the result
∞
q 2l+1
Vps (M) = P (cosψ )
q,M0 2 l(l +1) l M0M
l=1
X
q ψ
= qlogsin M0M , (19)
−2 − 2
from which the electric field is readily obtained
1
Eps (M) = Vps (M)
q,M0 −R∇S2 q,M0
q ψ
= cot M0M t (M) , (20)
2R 2 M0M
an expression which can also be obtained by applying Gauss theorem [4, 5]. Some remarks
are in order.
(i) One checks that, of course one has : Vbi (M) = Vps (M) Vps (M), since the
• q,M0 q,M0 − q,M0
backgrounds of the two pseudo-charges cancell each other.
(i) For R = 1 one has
• 6
1
∆ Vps (M) = 2πq δ(M ,M) . (21)
S2(O,R) q,M0 − { 0 − S}
where S = 4πR2 is the surface of (0,R).
2
S
(iii) Since the 3D distance between points M and M, i.e. the length of the chord
0
•
joining the two points is d = 2Rsin(ψ /2), it turns out that the expression of the
M0M
potential of a pseudo-charge in coincides with that of a point charge in the plane
2
S
E .
2
(iv) Although the potential of a pseudo-charge does not satisfies to Laplace-Beltrami
•
equation at source-free points, the potentials created by neutral multipoles made of
pseudo-charges indeed do, since the backgrounds cancell, as it will be seen below.
In order to extend the Green’s function expansion (4) to the sphere we have found the
following trick. Let us introduce, as in Ref. [27], the Cayley-Klein parameters
θ
α = exp(iϕ/2)cos , (22a)
2
θ
β = iexp( iϕ/2)sin . (22b)
− − 2
9
A short calculation reveals that, for two unit vectors z and z of , one has
1 2 2
S
ψ
12
α β α β = sin , (23)
1 2 2 1
| − | 2
from which it follows that
ψ θ θ
12 > <
logsin = logsin +logcos +log 1 zexp(i∆ϕ) , (24)
2 2 2 | − |
where θ = sup(θ ,θ ), θ = inf(θ ,θ ), ∆ϕ = (ϕ ϕ ), and z = tan(θ /2)/tan(θ /2).
> 1 2 < 1 2 1 2 < >
± −
We note that log 1 zexp(i∆ϕ) = (1/2)[log(1 zexp(i∆ϕ))+log(1 zexp( i∆ϕ))] and
| − | − − −
that, since z < 1, we can use the expansion of the complex logarithm log(1 Z) within its
| | −
circle of convergence, i.e. log(1 Z) = ∞ Zn/n for Z < 1, where Z = zexp( i∆ϕ).
− − n=1 | | ±
In doing so, we obtain the following expanPsion :
ψ θ θ
12 > <
logsin = logsin logcos
− 2 − 2 − 2
∞ n
1 tan(θ /2)
<
+ cos(n∆ϕ) , (25)
n tan(θ /2)
n=1 (cid:20) > (cid:21)
X
which seems to be a new mathematical result. In order to make something usefull of that
esthetic formula, we consider now a set of N pseudo-charges q of (O,R), with spherical
i 2
S
coordinates (θ ,ϕ ), and all around the north pole , i.e. π θ > θ , i = 1,...,N. A
i i 0 i
N ≥ ∀
multipolar expansion of the electrostatic potential Vps (M) created by these charges at
S2(O,R
some point M = (θ,ϕ) of the surface of the sphere away from follows from Eq. (25) under
N
the assumption that π θ θ . One finds
0
≥ ≥
N N
q ψ
Vps (M) = i q logsin MiM
S2(O,R) − 2 − i 2
i=1 i=1
X X
∞
1 θ 1 1
= Vps Q log 2Rsin + Q exp( inνϕ) , (26)
0 − 2 0ν 2 2 nXn nν −
ν (cid:20) (cid:21) n=1 ν
X XX
where Vps is some unessential constant. In (26) the index ν assumes, as in Sec. IIA, the
0
values 1 and +1 only. On the sphere, the complex multipolar moments are defined as
−
N
Q = q Xnexp(inνϕ ) . (27)
nν i i i
i=1
X
In Eqs. (26) and (27) we have introduced the variables X = 2Rtan(θ/2) not to be confused
⌢
with the length 2Rsin(θ/2) of the chord NM which constitutes the argument of the “log”
term in the r.h.s. of (26).
10
