Shape-Adapted Ewald Summation

M. Widom Department of Physics, Carnegie Mellon University, Pittsburgh, PA 15213

Ewald summation for electro- and magnetostatic interactions with replaces the conditionally convergent sum over distant charges with an absolutely convergent sum whose exact form depends on the choice of boundary conditions at infinity. We consider the Ewald summation in systems with periodic boundary conditions extending to infinity subject to a definite global sample shape. Our derivation, which explicitly considers ellipsoidal shapes, shows that a shape-dependent global depolarization energy supplements the pairwise interactions of particles. We derive explicit formulas for charge-charge and dipole-dipole interactions. The method is applied to Monte Carlo simulations of dipolar hard spheres to investigate shape-dependent magnetic susceptibility. It should be especially useful for systems in applied external fields and for systems that spontaneously polarize.

I. INTRODUCTION

Long-range electro- and magneto-static interactions present unusual challenges for theoretical analysis and computer simulation. In special circumstances (e.g. in the presence of applied fields or in the absence of charge neutrality) these systems lack conventional thermodynamic limits. Their free energies exhibit global system shape-dependence such as depolarization energies, and shape-dependent susceptibilities [1–4]. Many interesting ferrohydrodynamic behaviors of ferrofluids [5] arise from such effects.

Even cases with proper thermodynamic limits (e.g. charge neutrality, no applied fields [6–9]) require special handling of the long-range interaction to ensure correct results. For example, in the event of spontaneous polarization care is needed to remove depolarization energies from calculations. In Nature, this occurs spontaneously through the mechanism of domain formation [10]. Experimental studies achieve single-domain systems through the use of highly elongated needle-shaped samples. Even in analytical calculations must take proper account of this effect [11–13]

The following analysis concentrates on the case of Coulombic interactions between point charges. Other long-range interactions (specifically, charge-dipole, dipole-dipole, and charge-quadrupole) may be derived by taking limits of charge distributions, resulting in gradients of the Coulomb interaction.

We focus on the problem of computer simulations that utilize a small number of particles in an attempt to mimic the thermodynamic behavior of macroscopic systems. To avoid errors caused by large surface/volume ratios, simulations usually employ periodic boundary conditions. In a cubic box of edge length L, the interaction Hamiltonian including periodic boundaries in the repeated image convention is

H = \frac{1}{2} ∑_{i=1}^{N} ∑_{j=1}^{N} ∑_{n}^{′} \frac{q_{i}q_{j}}{|r + Ln|}  (1)

where r_{ij} = r_{i} −r_{j} is the separation of the particles i and j in the minimum image convention. The innermost sum extends over all simple cubic lattice points with integer coordinates, n = (n_{x}, n_{y}, n_{z}), and the prime on the sum indicates that n = 0 should be omitted in the case that i = j.

Owing to the periodic boundary conditions, there is no surface separating the simulated particles from a surrounding environment. However, if the interaction ϕ is long-ranged, falling off in d dimensions as 1/r^{d} or slower, the summations in (1) are conditionally convergent, so we must specify the manner in which the sum over n is to be performed. In particular, we must specify the outer shape of the set {n}, and the nature of the surrounding environment, even though these lie at infinity.

We consider the case in which the periodic boundaries are repeated throughout an infinitely large ellipsoidal shape and surrounded by vacuum. Our chief result rewrites the Hamiltonian in (1) as

H = \frac{1}{2} ∑_{i≠j=1}^{N} \frac{1}{L}q_{i}q_{j}ψ(\frac{r_{ij}}{L}) + ξ ∑_{i=1}^{N} q_{i}^{2} + \frac{2π}{L^{3}} P · D · P. (2)

which includes the usual result of Ewald summation,

ψ(r) = ∑_{n} \frac{erfc(α(|r + n|)}{|r + n|}+\frac{1}{π} ∑_{n≠0} |n|^{−2} exp [2πin · r − π^{2}|n|^{2}/α^{2}],

(3)

as well as the self-energy of a lattice of charges with

ξ = lim_{|r|→0}[ψ(r) − |r|^{−1}],  (4)

and finally, the surface dipole term in which

P = ∑_{i=1}^{N} q_{i}r_{i}  (5)

is the total system dipole moment, and D is the shapedependent depolarization tensor. This final term equals the depolarization energy arising from a surface pole density σ = P · n̂, where internal polarization at the surface has a component parallel to the surface normal n̂.

The depolarization tensor D takes special forms depending on the ellipsoid shape, such as

D_{s} = ( \begin{matrix}\frac{1}{3} & 0 & 0 \\ 0 & \frac{1}{3} & 0 \\ 0 & 0 & \frac{1}{3} \end{matrix} ) ,  (6)

for a sphere,

D_{n} = ( \begin{matrix}\frac{1}{2} & 0 & 0 \\ 0 & \frac{1}{2} & 0 \\ 0 & 0 & 0 \end{matrix} ) ,  (7)

for a needle-shape highly elongated in the z direction, and

D_{p} = ( \begin{matrix}0 & 0 & 0 \\ 0 & 0 & 0 \\ 0 & 0 & 1 \end{matrix} ) ,  (8)

respectively, for a sphere, a needle-shape highly elongated in the z direction, and for a pancake-shape flattened in the z direction

Our result differs from the classic result of de Leeuw, Perram and Smith [14] only in this final surface dipole term, which reduces to their expression in the special case of a spherical outer boundary shape. Typical simulations using Ewald summation drop the surface dipole term, because it limits the susceptibility and inhibits spontaneous polarization, obscuring the intrinsic thermodynamic properties. Dropping this term creates a “tin-foil” boundary condition, which wraps the surface at infinity in a material that neutralizes the surface pole density σ. Our result offers an alternative, in which a needle-like shape is taken so that one component of the depolarization tensor vanishes, eliminating depolarization fields in this direction, while maintaining them in the perpendicular directions. This choice is consistent with common experimental practice, where needle-shaped samples are preferred for measurements of susceptibility and polarization.

The following section derives eq. (2) for charge-charge interactions. All electrostatic and magnetostatic interactions follow from this basic result, as we show explicitly for dipole-dipole interactions in section 37. Finally, in section IV, we apply our result to simulate the shapedependent susceptibility of a dipolar fluid.

II. DERIVATION

The quadratic equation

r · E · r = \frac{1}{R^{2}} ,  (9)

with

E = ( \begin{matrix}\frac{1}{a^{2}} & 0 & 0 \\ 0 & \frac{1}{b^{2}} & 0 \\ 0 & 0 & \frac{1}{c^{2}} \end{matrix} )  (10)

defines an ellipsoid of semi-axes aR, bR, cR. We wish to evaluate the infinite lattice sum over n in the Hamiltonian (1) by summing first over the set of n values contained inside the ellipsoid of semi-axes a, b, c, then adding the ellipsoidal shell extending to semi-axes 2a, 2b, 2c, then extending to 3a, 3b, 3c, etc. Multiplying each term in the Hamiltonian (1) by a convergence factor

f(n, s) = e^{−sn·E·n}  (11)

achieves the same effect. When s > 0, the sum over ellipsoids becomes absolutely convergent, and the sum in (1) corresponds to smoothly cutting off the sum at the ellipsoid of semi-axes aR, bR, cR with R = 1/\sqrt{s}. The limit of s → 0 corresponds to extending the sum to infinity within this shape.

For charge-charge interactions the energy including the convergence factor becomes

E(s) = \frac{1}{2} ∑_{i=1}^{N} ∑_{j=1}^{N} q_{i}q_{j} ∑_{n}^{′} |r + n|^{−1}e^{−sn·E·n} (12)

= \frac{1}{2} ∑_{i≠j=1}^{N} q_{i}q_{j}ψ̄(r, s) + \frac{1}{2} ∑_{i=1}^{N} q_{i}^{2}ξ̄  (13)

where, again, the prime indicates dropping the n = 0 term when i = j. In the second line this interaction is broken into the effective interaction between pairs of charges at separation r (and all their periodic images),

ψ̄(r, s) = ∑_{n} \frac{e^{−sn·E·n}}{|r + n|}  (14)

and the effective interaction of individual charges with their own images,

ξ̄(s) = ∑_{n≠0} \frac{e^{−sn·E·n}}{|n|} .  (15)

The main steps of our calculation follow the method of de Leeuw, Perram and Smith [14] whose convergence factor f(n, s) = e^{−s|n|^{2}} corresponds to the special case of spherical shape. We omit many steps that are identical to theirs. First, replace |r+n|^{−1} in (14) using an integral representation of the Gamma-function Γ(\frac{1}{2}) to obtain

ψ̄(r, s) = \frac{1}{\sqrt{π}} ∑_{n} ∫_{0}^{∞} t^{− \frac{1}{2}} e^{−sn·E·n−t|r+n|^{2}}dt. (16)

Now break up the integral over t into a low t region 0 ≤ t ≤ α^{2} and a high t region α^{2} ≤ t ≤ ∞, so that

ψ̄(r, s) = ψ̄_{<}(r, s) + ψ̄_{>}(r, s)  (17)

The high t integral

ψ̄_{>}(r, s) ∑_{n} ∫_{α^{2}}^{∞} t^{− \frac{1}{2}} e^{−sn·E·n−t|r+n|^{2}}dt  (18)

is nonsingular and reduces to a sum of complementary error functions in the s → 0 limit.

The low t integral contains the singularity as s → 0 and requires additional work. Introduce matrix notation for the exponent in ψ̄,

F (r, n) ≡ sn · E · n + t(r + n) · I · (r + n)  (19)

with I the unit matrix. Defining U = sE + tI, rewrite

F (r, n) = G(r, n) + H(r)  (20)

with

G(r, n) = (n + tr · U^{−1}) · U · (n + tU^{−1} · r) (21)
H(r) = str · (U^{−1}E) · r

so that

ψ̄_{<}(r, s) = \frac{1}{\sqrt{π}} ∫_{0}^{α^{2}} t^{− \frac{1}{2}} e^{−H(r)} ∑_{n} e^{−G(r,n)}dt. (22)

The Jacobi Imaginary Transformation transforms the sum over n into

∑_{n} e^{−G(r,n)} = (\frac{π}{D})^{3/2} ∑_{n} e^{−π^{2}n·U^{−1}·n+2πitn·U^{−1}·r} (23)

with D^{3} = (t + s/a^{2})(t + s/b^{2})(t + s/c^{2}) the determinant of U.

The advantage of this new expression is that its n = 0 term contains the entire singularity of (23) as s → 0. Separate this term from the others, so that ψ̄_{<}(r, s) = ψ̄_{<}^{0}(r, s) + ψ̄_{<}^{≠0}(r, s) and note that

ψ(r) ≡ lim_{s→0}(ψ̄_{>}(r, s) + ψ̄_{<}^{≠0}(r, s))  (24)

takes the form given in (3).

It remains to evaluate

ψ̄_{<}^{0}(r, s) = π ∫_{0}^{α^{2}} \frac{e^{−H(r)}dt}{\sqrt{t(t + s/a^{2})(t + s/b^{2})(t + s/c^{2})}}

(25)

in the limit of small s. To do so, expand the exponential in powers of H, change variables from t to u = t/s, then integrate. Only the first two terms survive the s → 0 limit. In the first term,

ψ̄_{0}(r, s) = \frac{1}{s}ψ̄_{<}^{0,0}(s) + ψ̄_{<}^{0,1}(r, s),  (26)

where

lim_{s→0} ψ̄_{0,0}(s) = ∫_{0}^{∞} \frac{dt}{\sqrt{t(t + a^{2})(t + b^{2})(t + c^{2})}} ≡ \frac{abc}{C}

(27)

is inverse to the ellipsoid capacitance C [1]. In the second term,

lim_{s→0} ψ̄_{0,1}(r, s) ≡ −2πr · D · r  (28)

where

D = \frac{abc}{2} ∫_{0}^{∞} \frac{(tI + E)^{−1} · Edt}{\sqrt{t(t + a^{2})(t + b^{2})(t + c^{2})}} (29)

is the depolarization factor [1] of the ellipsoid.

Similar calculations for the self-interaction yield ξ̄ = ξ + \frac{ψ̄_{<}^{0,0}}{s}, with ξ as given in (4). Now collect the mutual and self interactions, insert into the energy equation (12),

and take the s → 0 limit. The terms proportional to \frac{ψ̄_{<}^{0,0}}{s} cancel owing to charge neutrality. Charge neutrality also yields

− ∑_{i≠j=1}^{N} q_{i}q_{j}r_{ij} · D · r_{ij} = P · D · P  (30)

Finally, collecting all terms and scaling lengths by L yields the result in eq. (2). Although the first two terms in (2) depend on the parameter α, the value of the Hamiltonian H does not.

III. DIPOLE-DIPOLE INTERACTIONS

Systems of particles interacting through dipole-dipole interactions provide a fascinating range of interesting behavior, including chain formation at low densities [15– 18] to possible spontaneous polarization in the liquid state [19, 20]. Spontaneous polarization is observed when the surface dipole term is removed using tin-foil boundaries, otherwise polarized domains arise [21, 22]

We model the dipole-dipole interaction as a limit of the charge-charge interaction in which each charge q_{i} at r_{i} (i = 1, …, N) pairs with an opposite charge q_{i+N} at an infinitesimally close position r_{i+N} = r_{i} + δ_{i}, creating a point dipole moment μ_{i} = −q_{i}δ_{i}. The first term in the charge-charge interaction (2) breaks into four distinct contributions

\frac{1}{2} ∑_{i≠j=1}^{2N} q_{i}q_{j}ψ(r_{ij}) = \frac{1}{2} ∑_{i≠j=1}^{N} q_{i}q_{j}ψ(r_{ij}) + \frac{1}{2} ∑_{i≠j=1}^{N} q_{i}q_{j}ψ(r_{ij} + δ_{i} − δ_{j}) (31)

−\frac{1}{2} ∑_{i=1}^{N} q_{i}^{2}ψ(δ_{i}) − \frac{1}{2} ∑_{i≠j=1}^{N} q_{i}q_{j}ψ(r_{ij} + δ_{i}).

To this expression we must add the charge lattice selfenergy ξ ∑_{i=1}^{N} q_{i}^{2} and the surface dipole term 2πP·D·P, with P = ∑_{i=1}^{N} μ_{i}.

We approximate terms in the sum (31 above using expansions for small δ. First, expand

ψ(δ) = \frac{1}{δ} + ξ + \frac{2π}{3} δ^{2} + ⋯  (32)

The leading term in this expansion corresponds to the direct interaction between charges q_{i}q_{i+N}/δ forming each dipole μ_{i}. This energy should be allocated to the formation energy for isolated dipoles, and hence we remove it from the total energy. The next order in the expansion contributes q_{i}^{2}ξ. These terms cancel against the charge lattice self-energy. The final term in the expansion contributes

−\frac{2π}{3} |μ_{i}|^{2}  (33)

to the total energy, which we identify as the dipole lattice self-energy. Next, expand

ψ(r+δ) = ψ(r)+(δ ·∇)ψ(r)+ \frac{1}{2}(δ ·∇)^{2}ψ(r)+· · · . (34)

The resulting terms of form q_{i}q_{j}ψ(r_{ij}) cancel each other in (31), as do also the terms of form q_{i}q_{j}(δ · ∇)ψ(r). The quadratic term in this expansion contributes

− ∑_{i≠j=1}^{N}(μ_{i} · ∇)(μ_{j} · ∇)ψ(r)  (35)

to the total energy, which we identify as the dipole-dipole interaction energy.

We have now fully accounted for the first two terms in the charge-charge Hamiltonian (2). The final term, the surface dipole, separates into a surface self-interaction 2πμ_{i} · D · μ_{i}, and a mutual interaction 2πμ_{i} · D · μ_{j}. Collecting all terms, our final expression for the total energy of a collection of dipole moments is

H = H_{self} + H_{pair}  (36)

where

H_{self} = −2π ∑_{i=1}^{N} μ_{i} · (D − \frac{1}{3}I) · μ_{i}  (37)

H_{pair} = \frac{1}{2} ∑_{i≠j=1}^{N}[4πμ_{i} · D · μ_{j} − (μ_{i} · ∇)(μ_{j} · ∇)ψ(r_{ij})].

The classic result of de Leeuw, Perram and Smith [14] deals with the case of spherical sample, for which D = \frac{1}{3}I and thus H_{self} vanishes.

IV. SHAPE-DEPENDENT SUSCEPTIBILITY OF DIPOLAR FLUID

As a first application of our new Hamiltonian, we simulate the magnetic susceptibility of a dipolar fluid assuming a highly elongated needle-like sample shape, with depolarization factor D_{n} given in (7). In a general ellipsoidal sample, the observed magnetic susceptibility is [3]

χ_{obs} = \frac{χ_{int}}{1 + 4πDχ_{int}}  (38)

where χ_{int} is an intrinsic susceptibility representing the polarization response to an applied field in the absence of surface depolarization effects, and D is the appropriate element of the depolarization factor. In the special case of a needle-like shape lying along the z axis,

χ_{xy} = \frac{χ_{int}}{1 + 2πχ_{int}}  χ_{z} = χ_{int}  (39)

We choose the dipolar hard sphere model as the simplest representation of a dipolar fluid. Such models have been employed to crudely describe electrically polar fluids such as H_{2}O, and magnetic fluids such as ferrofluid. Introducing the hard sphere diameter a, a dimensionless measure of dipolar coupling strength is λ = μ^{2}/k_{B}T a^{3} and we define a dimensionless temperature T^{⋆} = 1/λ. For a simulation of N particles in a box of volume V we introduce a dimensionless density ρ^{⋆} = Na^{3}/V. The thermodynamic properties of the fluid depend only on T^{⋆} and ρ^{⋆}.

The limit of very low density obeys a Virial expansion [15, 23–25]. Here we consider the expansions of pressure and susceptibility in powers of ρ^{⋆},

\frac{1}{k_{B}T} P (T^{⋆}, ρ^{⋆}) = B_{1}ρ^{⋆} + B_{2}ρ^{⋆}^{2} + B_{3}ρ^{⋆}^{3} + ⋯ (40)
χ_{int}(T^{⋆}, ρ^{⋆}) = M_{1}ρ^{⋆} + M_{2}ρ^{⋆}^{2} + M_{3}ρ^{⋆}^{3} + · · · .

Thermodynamic identities relate the susceptibility coefficients M_{n} to magnetic field derivatives of the pressure coefficients B_{n}. In the limit of low density, the pressure and susceptibility takes the classical ideal gas values

B_{1} = 1  M_{1} = \frac{1}{3T^{⋆}} .  (41)

The higher order terms have simple representations in the limits of strong or weak coupling. We reproduce here the strong coupling limits, because the third-order terms [24] have not previously been published:

\begin{matrix}B_{2} = −\frac{πa^{3}e^{2λ}}{18λ^{3}} & M_{2} = −\frac{2}{3}λa^{3}B_{2} \\ B_{3} = −\frac{32768π^{2}a^{6}e^{17λ/4}}{7373835λ^{6}} & M_{3} = −λa^{3}B_{3}. \end{matrix}  (42)

In the strong coupling limit, the virial coefficients reflect the association of dipolar particles into chains [16]. For example, the densities of two- and three-particle chains are ρ_{2} = −B_{2}ρ^{⋆} and ρ_{3} = −\frac{1}{2}B_{3}ρ^{⋆}^{2}, while the number of isolated particles ρ_{1} = ρ^{⋆} − 2ρ_{2} − 3ρ_{3}. The susceptibility coefficients then follow from treating the length-n chains as effective dipoles of strength nμ. Namely [24, 26],

χ = \frac{1}{3T^{⋆}} ρ_{1} + \frac{2^{2}}{3T^{⋆}} ρ_{2} + \frac{3^{2}}{3T^{⋆}} ρ_{3} + · · · .  (43)

The interesting shape-dependent effects enter as χ_{z} grows large. Hence we focus our simulation on relatively low temperature T^{⋆} = 0.2 and follow a sequence of increasing density. At this low temperature particles tend to associate into chains clusters of varying size. Fig. 1 illustrates a typical configuration. Simulations were performed for systems of N = 100 particles in cubic boxes of edge length L = 10 − 100. Chain formation causes slow evolution of moment orientations through single particle rotations, reducing the efficiency with which we sample configuration space. To counteract this effect, we employ moment reversal within entire chains as one of our Monte Carlo steps [27].

Fig. 2 illustrates our data for χ_{xy} and χ_{z}. The shape effect is clearly observed, as χ_{z} grows without bound while χ_{xy} saturates at 1/2π = 0.159. The first, second and third order Virial results are included for comparison, showing excellent agreement at low densities but growing inaccurate as particles associate into long chains and branched structures. In particular, the virial series tends to overestimate χ_{z} because up to third order it only includes the high susceptibility contributions of linear chains, while low susceptibility closed rings and other structures become favorable for lengths of 5 or more particles [28].

In conclusion, we have derived the Ewald summation formula for electro- and magneto-static interactions that utilize periodic boundary conditions in a cubic box, but where the periodic replicas extend to infinity within a global ellipsoidal shape. The resulting shape-dependent surface dipole term creates depolarization effects. In the limit of needle-like shapes the depolarization fields vanish parallel to the long axis. Such a sample shape should be ideal for simulations in the presence of applied fields or spontaneous polarization, and are more consistent with actual experimental geometries than alternative approaches such as tin-foil boundary conditions.

V. ACKNOWLEDGMENTS

This work was supported in part by NSF grant DMR-0111198. Useful discussions with R. B. Griffiths are acknowledged.

[1] L. D. Landau and E. M. Lifshitz, Electrodynamics of continuous media (Pergamon Press, Oxford, 1960)

[2] W. F. Brown, Magnetostatic principles in ferromagnetism (North-Holland, Amsterdam, 1962)

[3] J. A. Osborn, Phys. Rev. 67, 351 (1945)

[4] J. M. Luttinger and L. Tisza, Phys. Rev. 70 954 (1946)

[5] R. E. Rosensweig, Ferrohydrodynamics (Cambridge, 1985)

[6] R. B. Griffiths, Phys. Rev. 176, 655 (1968)

[7] E. H. Lieb, Rev. Mod. Phys. 48, 553 (1976)

[8] S. Banerjee, R. B. Griffiths and M. Widom, J. Stat. Phys. 93, 109 (1998)

[9] S. Banerjee, R. B. Griffiths and M. Widom, J. Stat. Phys. 104, 725 (1998)

[10] F. Bitter, Phys. Rev. 38, 1903 (1931)

[11] H. Zhang and M. Widom, Phys. Rev. E49, R3591 (1994)

[12] B. Groh and S. Dietrich, Phys. Rev. Lett. 72, 2422 (1994)

[13] M. Widom and H. Zhang, Phys. Rev. Lett. 74, 2616 (1995)

[14] S. W. de Leeuw, J. W. Perram and E. R. Smith, Proc. Roy. Soc. Lond. A 373, 27 (1980)

[15] P. G. de Gennes and P. A. Pincus, Phys. Kond. Mater. 11, 189 (1970)

[16] P. C. Jordan, Mol. Phys. 25, 961 (1973)

[17] J. J. Weis and D. Levesque, Phys. Rev. E 49, 5131 (1993)

[18] M. E. Van Leeuwen and B. Smit, Phys. Rev. Lett.71, 3991 (1993)

[19] D. Wei and G. Patey, Phys. Rev. Lett.68, 2043 (1992)

[20] J. J. Weis, D. Levesque and G. J. Zarragoicoechea, Phys. Rev. Lett. 69, 913 (1992)

[21] D. Wei and G. Patey, Phys. Rev. A46, 7783 (1992)

[22] J. J. Weis and D. Levesque, Phys. Rev. E48, 3728 (1992)

[23] C. Joslin and S. Goldman, Mol. Phys. 79, 499 (1993)

[24] H. Zhang, thesis (Carnegie Mellon University, 1994)

[25] B. Huke and M. Lucke, Phys. Rev. E62, 6875 (2000)

[26] M. Widom and H. Zhang, Mat. Res. Soc. Symp. 248, 235 (1992)

[27] S. W. Davis, W. McCausland, H. C. McGahagan, C. T. Tanaka and M. Widom, Phys. Rev. E59, 2424 (1999)

[28] P. Jund, S. G. Kim, D. Tomanek and J. Hetherington, Phys. Rev. Lett. 74, 3049 (1995)

FIG. 1: (left) Typical configuration at T^{⋆} = 0.2, ρ^{⋆} = 0.0125 projected into yz plane. Black and white dots indicate north and south magnetic poles.

FIG. 2: Susceptibility data versus ρ^{⋆} at fixed T^{⋆} = 0.2. Top graph illustrates intrinsic susceptibility χ_{int} = χ_{z}, bottom graph illustrates χ_{xy}. Solid, dotted and dashed lines, respectively, are first, second and third order Virial series in the strong coupling limit.