\documentclass[twocolumn,amsmath,amssymb]{revtex4}
\usepackage{graphicx}% Include figure files
\usepackage{dcolumn}% Align table columns on decimal point
\usepackage{bm}% bold math
\def\be{\begin{equation}}
\def\ee{\end{equation}}
\def\bn{{\bf n}}
\def\br{{\bf r}}
\def\bD{{\bf D}}
\def\bI{{\bf I}}
\def\bN{{\bf N}}
\def\bP{{\bf P}}
\def\bU{{\bf U}}
\def\bE{{\bf E}}
\def\rs{{\rho^{\star}}}
\def\ts{{T^{\star}}}
\def\bmu{{\mbox{\boldmath $\mu$}}}
\def\bdelta{{\mbox{\boldmath $\delta$}}}
\begin{document}

%\preprint{cond-mat/04xxxxxx}

\title{Shape-Adapted Ewald Summation}

\author{M. Widom
\newline
Department of Physics, Carnegie Mellon University,
Pittsburgh, PA  15213
}

\begin{abstract}
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.
\end{abstract}

%\pacs{Valid PACS appear here}


\maketitle

\section{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~\cite{LL,Brown,Osborn45,LT46}. Many interesting
ferrohydrodynamic behaviors of ferrofluids~\cite{RR} arise from such
effects.

Even cases with proper thermodynamic limits (e.g. charge neutrality,
no applied fields~\cite{RBG68,Lieb76,BGW98,BGW01}) 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~\cite{Bitter31}.  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~\cite{ZW94,GD94,WZ95}

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
\begin{equation}
H = \frac{1}{2} \sum_{i=1}^N \sum_{j=1}^N \sum_{{\bn}}^{\prime}
\frac{q_i q_j}{|\br+L\bn|}
\label{Ham}
\end{equation}
where $\br_{ij}=\br_i-\br_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, $\bn=(n_x,n_y,n_z)$, and the prime on the sum indicates that
$\bn={\bf 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 $\phi$ is long-ranged, falling off in $d$
dimensions as $1/r^d$ or slower, the summations in (\ref{Ham}) are
conditionally convergent, so we must specify the manner in which the sum
over $\bn$ is to be performed.  In particular, we must specify the outer
shape of the set $\{\bn\}$, and the nature of the surrounding
environment, {\em 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 (\ref{Ham}) as
\begin{equation}
H = \frac{1}{2}\sum_{i\ne j=1}^N \frac{1}{L} q_iq_j \psi(\frac{\br_{ij}}{L})
+ \xi \sum_{i=1}^N q_i^2 
+ \frac{2\pi}{L^3} \bP\cdot\bD\cdot\bP.
\label{Ewald}
\end{equation}
which includes the usual result of Ewald summation,
\begin{equation}
\label{psi}
\psi(\br) =  \sum_{\bn} \frac{{\rm erfc}(\alpha(|\br+\bn|)}{|\br+\bn|}
 +
\frac{1}{\pi} \sum_{\bn\ne {\bf 0}}
|\bn|^{-2} \exp{[2\pi i \bn \cdot \br - \pi^2|\bn|^2/\alpha^2]},
\end{equation}
as well as the {\em self-energy} of a lattice of charges with
\begin{equation}
\xi = \lim_{|\br|\rightarrow 0} [\psi(\br) - |\br|^{-1}],
\label{xi}
\end{equation}
and finally, the {\em surface dipole} term in which
\begin{equation}
\bP = \sum_{i=1}^N q_i \br_i
\label{P}
\end{equation}
is the total system dipole moment, and $\bD$ is the shape-dependent
depolarization tensor.  This final term equals the depolarization
energy arising from a surface pole density $\sigma=\bP\cdot\hat{n}$,
where internal polarization at the surface has a component parallel to
the surface normal $\hat{n}$.

The depolarization tensor $\bD$ takes special forms depending on the
ellipsoid shape, such as
\begin{equation}
\label{sphere}
\bD_s = \left(
\begin{array}{ccc}
\frac{1}{3} & 0 & 0 \\
0 & \frac{1}{3} & 0 \\
0 & 0 & \frac{1}{3}
\end{array}
\right),
\end{equation}
for a sphere,
\begin{equation}
\label{needle}
\hspace{1cm}
\bD_n = \left(
\begin{array}{ccc}
\frac{1}{2} & 0 & 0 \\
0 & \frac{1}{2} & 0 \\
0 & 0 & 0
\end{array}
\right),
\end{equation}
for a needle-shape highly elongated in the $z$ direction, and
\begin{equation}
\label{pancake}
\hspace{1cm}
\bD_p = \left(
\begin{array}{ccc}
0 & 0 & 0 \\
0 & 0 & 0 \\
0 & 0 & 1
\end{array}
\right),
\end{equation}
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~\cite{dLPS80} 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 $\sigma$.  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.~(\ref{Ewald}) for charge-charge
interactions.  All electrostatic and magnetostatic interactions follow
from this basic result, as we show explicitly for
dipole-dipole interactions in section~\ref{dip-dip}.  Finally, in
section~\ref{Sim}, we apply our result to simulate the shape-dependent
susceptibility of a dipolar fluid.

\section{Derivation}
\label{Derive}

The quadratic equation
\begin{equation}
\br \cdot \bE \cdot \br = \frac{1}{R^2},
\end{equation}
with
\begin{equation}
\bE = \left(
\begin{array}{ccc}
\frac{1}{a^2} & 0 & 0 \\
0 & \frac{1}{b^2} & 0 \\
0 & 0 & \frac{1}{c^2}
\end{array}
\right)
\end{equation}
defines an ellipsoid of semi-axes $aR,bR,cR$.  We wish to
evaluate the infinite lattice sum over $\bn$ in the Hamiltonian
(\ref{Ham}) by summing first over the set of $\bn$ 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 (\ref{Ham}) by a convergence factor
\begin{equation}
f(\bn,s)=e^{-s \bn \cdot \bE \cdot \bn}
\end{equation}
achieves the same effect.  When $s>0$, the sum over ellipsoids becomes
absolutely convergent, and the sum in~(\ref{Ham}) 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\rightarrow 0$ corresponds to
extending the sum to infinity within this shape.

For charge-charge interactions the energy including the convergence
factor becomes
\begin{eqnarray}
\label{Eofs}
E(s) & = & \frac{1}{2} \sum_{i=1}^N \sum_{j=1}^N q_i q_j
\sum_{\bn}^{\prime} |\br+\bn|^{-1} e^{-s \bn \cdot \bE \cdot \bn} \\
& = & \frac{1}{2} \sum_{i\ne j=1}^N q_i q_j \bar{\psi}(\br,s)
+ \frac{1}{2} \sum_{i=1}^N q_i^2 \bar{\xi}
\end{eqnarray}
where, again, the prime indicates dropping the $\bn={\bf 0}$ term when
$i=j$.  In the second line this interaction is broken into the
effective interaction between pairs of charges at separation $\br$
(and all their periodic images),
\begin{equation}
\bar{\psi}(\br,s)=\sum_{\bn} \frac{e^{-s\bn\cdot\bE\cdot\bn}}{|\br+\bn|}
\label{psibar}
\end{equation}
and the effective interaction of individual charges with their
own images,
\begin{equation}
\bar{\xi}(s)=\sum_{\bn\ne{\bf 0}} \frac{ e^{-s\bn\cdot\bE\cdot\bn}}{|\bn|}.
\end{equation}

The main steps of our calculation follow the method of de Leeuw,
Perram and Smith~\cite{dLPS80} whose convergence factor
$f(\bn,s)=e^{-s|\bn|^2}$ corresponds to the special case of spherical
shape.  We omit many steps that are identical to theirs.  First,
replace $|\br+\bn|^{-1}$ in (\ref{psibar}) using an integral
representation of the Gamma-function $\Gamma(\frac{1}{2})$ to obtain
\begin{equation}
\bar{\psi}(\br,s) = \frac{1}{\sqrt{\pi}}
\sum_{\bn} \int_0^{\infty} t^{-\frac{1}{2}}
e^{-s \bn \cdot \bE \cdot \bn - t |\br+\bn|^2} {\rm d}t.
\end{equation}
Now break up the integral over $t$ into a low $t$ region $0\le t \le
\alpha^2$ and a high $t$ region $\alpha^2 \le t \le \infty$, so that
\begin{equation}
\bar{\psi}(\br,s)=\bar{\psi}_<(\br,s)+\bar{\psi}_>(\br,s)
\end{equation}
The high $t$ integral
\begin{equation}
\bar{\psi}_>(\br,s) \sum_{\bn} \int_{\alpha^2}^{\infty} t^{-\frac{1}{2}}
e^{-s \bn\cdot\bE\cdot\bn - t |\br+\bn|^2} {\rm d}t 
\end{equation}
is nonsingular and reduces to a sum of complementary error functions
in the $s\rightarrow 0$ limit.

The low $t$ integral contains the singularity as $s\rightarrow
0$ and requires additional work.  Introduce matrix notation for the
exponent in $\bar{\psi}$,
\begin{equation}
F(\br,\bn) 
\equiv s \bn \cdot \bE \cdot \bn + t (\br+\bn) \cdot \bI \cdot (\br+\bn)
\end{equation}
with $\bI$ the unit matrix.  Defining $\bU = s\bE + t\bI$, rewrite
\begin{equation}
F(\br,\bn) = G(\br,\bn) + H(\br)
\end{equation}
with
\begin{eqnarray}
G(\br,\bn) & = &
(\bn + t \br\cdot\bU^{-1})\cdot\bU\cdot(\bn + t \bU^{-1}\cdot\br) \\ \nonumber
H(\br) & = & s t \br\cdot(\bU^{-1}\bE)\cdot\br
\end{eqnarray}
so that
\begin{equation}
\bar{\psi}_<(\br,s) = \frac{1}{\sqrt{\pi}}
\int_0^{\alpha^2} t^{-\frac{1}{2}}e^{-H(\br)}
\sum_{\bn} e^{-G(\br,\bn)} {\rm d}t.
\end{equation}
The Jacobi Imaginary
Transformation transforms the sum over $\bn$ into
\begin{equation}
\label{JIT}
\sum_{\bn} e^{-G(\br,\bn)} = (\frac{\pi}{D})^{3/2}
\sum_{\bn} e^{-\pi^2 \bn \cdot \bU^{-1} \cdot \bn
+ 2\pi i t \bn \cdot\bU^{-1}\cdot\br}
\end{equation}
with $D^3=(t+s/a^2)(t+s/b^2)(t+s/c^2)$ the determinant of $\bU$.

The advantage of this new expression is that its $\bn = {\bf 0}$ term
contains the entire singularity of~(\ref{JIT}) as $s\rightarrow 0$.
Separate this term from the others, so that $\bar{\psi}_<(\br,s)
=\bar{\psi}_<^0(\br,s)+\bar{\psi}_<^{\ne 0}(\br,s)$ and note that
\begin{equation}
\psi(\br) \equiv \lim_{s\rightarrow 0}
(\bar{\psi}_>(\br,s)+\bar{\psi}_<^{\ne 0}(\br,s))
\end{equation}
takes the form given in (\ref{psi}).

It remains to evaluate
\begin{equation}
\bar{\psi}_<^0(\br,s) = \pi \int_0^{\alpha^2} \frac{e^{-H(\br)} dt}
{\sqrt{t(t+s/a^2)(t+s/b^2)(t+s/c^2)}}
\end{equation}
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\rightarrow 0$ limit.  In the first
term,
\begin{equation}
\bar{\psi}_0(\br,s) = \frac{1}{s} \bar{\psi}_<^{0,0}(s) +
\bar{\psi}_<^{0,1}(\br,s),
\end{equation}
where
\begin{equation}
\lim_{s\rightarrow 0} \bar{\psi}_{0,0}(s) = \int_0^{\infty}
\frac{dt}{\sqrt{t(t+a^2)(t+b^2)(t+c^2)}} \equiv \frac{abc}{C}
\end{equation}
is inverse to the ellipsoid capacitance $C$~\cite{LL}.  In the second term,
\begin{equation}
\lim_{s\rightarrow 0} \bar{\psi}_{0,1}(\br,s) \equiv -2\pi \br\cdot\bD\cdot\br
\end{equation}
where
\begin{equation}
\bD = \frac{abc}{2} \int_0^{\infty}\frac{(t\bI + \bE)^{-1}\cdot\bE {\rm d}t}
{\sqrt{t(t+a^2)(t+b^2)(t+c^2)}}
\end{equation}
is the depolarization factor~\cite{LL} of the ellipsoid.

Similar calculations for the self-interaction yield
$\bar{\xi}=\xi+\frac{\bar{\psi}_<^{0,0}}{s}$, with $\xi$ as given
in~(\ref{xi}).  Now collect the mutual and self interactions, insert
into the energy equation~(\ref{Eofs}), and take the $s\rightarrow 0$
limit.  The terms proportional to $\frac{\bar{\psi}_<^{0,0}}{s}$
cancel owing to charge neutrality.  Charge neutrality also yields
\begin{equation}
-\sum_{i\ne j=1}^N q_i q_j \br_{ij}\cdot\bD\cdot\br_{ij}
=\bP\cdot\bD\cdot\bP
\end{equation}
Finally, collecting all terms and scaling lengths by $L$ yields the
result in eq.~(\ref{Ewald}).  Although the first two terms in
(\ref{Ewald}) depend on the parameter $\alpha$, the value of the
Hamiltonian $H$ does not.

\section{Dipole-dipole interactions}
\label{Dip-dip}

Systems of particles interacting through dipole-dipole interactions
provide a fascinating range of interesting behavior, including chain
formation at low densities~\cite{dGP70,J73,WL94,vLS93} to possible
spontaneous polarization in the liquid state~\cite{WP92,WLZ92}.
Spontaneous polarization is observed when the surface dipole term is
removed using tin-foil boundaries, otherwise polarized domains
arise~\cite{WP92a,WL93}

We model the dipole-dipole interaction as a limit of the charge-charge
interaction in which each charge $q_i$ at $\br_i$ ($i=1,\dots,N)$
pairs with an opposite charge $q_{i+N}$ at an infinitesimally close
position $\br_{i+N}=\br_i+\bdelta_i$, creating a point dipole moment
$\bmu_i=-q_i\bdelta_i$.  The first term in the charge-charge
interaction~(\ref{Ewald}) breaks into four distinct contributions
\begin{eqnarray}
\label{DD1}
\frac{1}{2} \sum_{i\ne j=1}^{2N} q_i q_j \psi(\br_{ij}) & = &
\frac{1}{2} \sum_{i\ne j=1}^N q_i q_j \psi(\br_{ij})
+ \frac{1}{2} \sum_{i\ne j=1}^N q_i q_j \psi(\br_{ij}+\bdelta_i-\bdelta_j)
\\  \nonumber
& & -\frac{1}{2} \sum_{i=1}^N q_i^2 \psi(\bdelta_i)
-\frac{1}{2} \sum_{i\ne j=1}^N q_i q_j \psi(\br_{ij}+\bdelta_i).
\end{eqnarray}
To this expression we must add the charge lattice self-energy $\xi
\sum_{i=1}^N q_i^2$ and the surface dipole term $2\pi\bP\cdot\bD\cdot\bP$,
with $\bP=\sum_{i=1}^N \bmu_i$.

We approximate terms in the sum (\ref{DD1} above using expansions for
small $\bdelta$.  First, expand
\begin{equation}
\psi(\bdelta) = \frac{1}{\delta} + \xi + \frac{2\pi}{3} \delta^2 + \cdots
\end{equation}
The leading term in this expansion corresponds to the direct
interaction between charges $q_i q_{i+N}/\delta$ forming each dipole
$\bmu_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 \xi$.  These terms
cancel against the charge lattice self-energy.  The final term
in the expansion
contributes 
\begin{equation}
-\frac{2\pi}{3}|\bmu_i|^2
\end{equation}
to the total energy, which we identify as the {\em dipole lattice
self-energy}.  Next, expand
\begin{equation}
\psi(\br+\delta) = \psi(\br) + (\bdelta\cdot\nabla) \psi(\br)
+ \frac{1}{2} (\bdelta\cdot\nabla)^2 \psi(\br)+\cdots.
\end{equation}
The resulting terms of form $q_i q_j\psi(\br_{ij})$ cancel each other
in~(\ref{DD1}), as do also the terms of form $q_i q_j
(\bdelta\cdot\nabla)\psi(\br)$.  The quadratic term in this expansion
contributes
\begin{equation}
-\sum_{i\ne j=1}^N (\bmu_i\cdot\nabla)(\bmu_j\cdot\nabla)\psi(\br)
\end{equation}
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~(\ref{Ewald}).  The final term, the surface
dipole, separates into a surface self-interaction
$2\pi\bmu_i\cdot\bD\cdot\bmu_i$, and a mutual interaction
$2\pi\bmu_i\cdot\bD\cdot\bmu_j$.  Collecting all terms, our final
expression for the total energy of a collection of dipole moments is
\begin{equation}
H = H_{self}+H_{pair}
\end{equation}
where
\begin{eqnarray}
\label{dip-dip}
H_{self} & = & -2\pi\sum_{i=1}^N \bmu_i\cdot(\bD-\frac{1}{3}\bI)\cdot\bmu_i \\
H_{pair} & = & \frac{1}{2} \sum_{i\ne j=1}^N
[4\pi\bmu_i\cdot\bD\cdot\bmu_j
-(\bmu_i\cdot\nabla)(\bmu_j\cdot\nabla)\psi(\br_{ij})]\nonumber.
\end{eqnarray}
The classic result of de Leeuw, Perram and Smith~\cite{dLPS80} deals with
the case of spherical sample, for which $\bD=\frac{1}{3}\bI$ and thus
$H_{self}$ vanishes.

\section{Shape-dependent susceptibility of dipolar fluid}
\label{Sim}

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 $\bD_n$ given
in~(\ref{needle}).  In a general ellipsoidal sample, the observed
magnetic susceptibility is~\cite{Osborn45}
\begin{equation}
\chi_{obs}=\frac{\chi_{int}}{1+4\pi D\chi_{int}}
\end{equation}
where $\chi_{int}$ is an {\em 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,
\begin{equation}
\label{Demag}
\chi_{xy}=\frac{\chi_{int}}{1+2\pi\chi_{int}} \hspace{1cm}
\chi_z=\chi_{int}
\end{equation}

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
$\lambda=\mu^2/k_BTa^3$ and we define a dimensionless temperature
$\ts=1/\lambda$.  For a simulation of $N$ particles in a box of volume
$V$ we introduce a dimensionless density $\rs=N a^3/V$.  The
thermodynamic properties of the fluid depend only on $\ts$ and $\rs$.

The limit of very low density obeys a Virial
expansion~\cite{dGP70,Joslin93,Zhang94,Huke00}.  Here we consider the
expansions of pressure and susceptibility in powers of $\rho^{\star}$,
\begin{eqnarray}
\label{Virial}
\frac{1}{k_BT}P(\ts,\rs)&=&B_1 \rs+B_2\rs^2+B_3\rs^3+\cdots \\
\chi_{int}(\ts,\rs)&=&M_1 \rs+M_2\rs^2+M_3\rs^3+\cdots. \nonumber
\end{eqnarray}
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
\begin{equation}
B_1=1 \hspace{1cm} M_1=\frac{1}{3\ts}.
\end{equation}
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~\cite{Zhang94} have not
previously been published:
\begin{equation}
\begin{array}{lcl}
B_2=-\frac{\pi a^3 e^{2\lambda}}{18\lambda^3} &\hspace{1cm}&
M_2=-\frac{2}{3}\lambda a^3 B_2 \\
B_3=-\frac{32768 \pi^2 a^6 e^{17\lambda/4}}{7373835 \lambda^6} &\hspace{1cm}&
M_3 = -\lambda a^3 B_3.
\end{array}
\end{equation}
In the strong coupling limit, the virial coefficients reflect the
association of dipolar particles into chains~\cite{J73}.  For
example, the densities of two- and three-particle chains are
$\rho_2=-B_2\rs$ and $\rho_3=-\frac{1}{2}B_3 \rs^2$, while the number
of isolated particles $\rho_1=\rs-2\rho_2-3\rho_3$.  The
susceptibility coefficients then follow from treating the length-$n$
chains as effective dipoles of strength $n\mu$.  Namely~\cite{WZ92,Zhang94},
\begin{equation}
\chi=\frac{1}{3\ts}\rho_1+\frac{2^2}{3\ts}\rho_2+
\frac{3^2}{3\ts}\rho_3+\cdots.
\end{equation}

The interesting shape-dependent effects enter as $\chi_z$ grows large.
Hence we focus our simulation on relatively low temperature $\ts=0.2$
and follow a sequence of increasing density.  At this low temperature
particles tend to associate into chains clusters of varying size.
Fig.~\ref{config} 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~\cite{Davis99}.

Fig.~\ref{data} illustrates our data for $\chi_{xy}$ and $\chi_z$.
The shape effect is clearly observed, as $\chi_z$ grows without bound
while $\chi_{xy}$ saturates at $1/2\pi=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 $\chi_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~\cite{Tomanek95}.



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.

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

%\bibliography{dip}
\begin{thebibliography}{9}

\bibitem{LL}
L. D. Landau and E. M. Lifshitz, {\em Electrodynamics of continuous media}
(Pergamon Press, Oxford, 1960)

\bibitem{Brown}
W. F. Brown, {\em Magnetostatic principles in ferromagnetism}
(North-Holland, Amsterdam, 1962)

\bibitem{Osborn45}
J. A. Osborn, {\em Phys. Rev.} {67}, 351 (1945)

\bibitem{LT46}
J. M. Luttinger and L. Tisza, {\em Phys. Rev.} {\bf 70} 954 (1946)

\bibitem{RR}
R. E. Rosensweig, {\em Ferrohydrodynamics} (Cambridge, 1985)

\bibitem{RBG68}
R. B. Griffiths, {\em Phys. Rev.} {\bf 176}, 655 (1968)

\bibitem{Lieb76}
E. H. Lieb, {\em Rev. Mod. Phys.} {48}, 553 (1976)

\bibitem{BGW98}
S. Banerjee, R. B. Griffiths and M. Widom, {\em J. Stat. Phys.} {\bf 93},
109 (1998)

\bibitem{BGW01}
S. Banerjee, R. B. Griffiths and M. Widom, {\em J. Stat. Phys.} {\bf 104},
725 (1998)

\bibitem{Bitter31}
F. Bitter, {\em Phys. Rev.} {\bf 38}, 1903 (1931)

\bibitem{ZW94}
H. Zhang and M. Widom, {\em Phys. Rev. E}{\bf 49}, R3591 (1994)

\bibitem{GD94}
B. Groh and S. Dietrich, {\em Phys. Rev. Lett.} {\bf 72}, 2422 (1994)

\bibitem{WZ95}
M. Widom and H. Zhang, {\em Phys. Rev. Lett.} {\bf 74}, 2616 (1995)

\bibitem{dLPS80}
S. W. de Leeuw, J. W. Perram and E. R. Smith, {\em Proc. Roy. Soc. Lond.}
{\bf A 373}, 27 (1980)

\bibitem{dGP70}
P. G. de Gennes and P. A. Pincus, {\em Phys. Kond. Mater.} {\bf 11}, 189 (1970)

\bibitem{J73}
P. C. Jordan, {\em Mol. Phys.} {\bf 25}, 961 (1973)

\bibitem{WL94}
J. J. Weis and D. Levesque, {\em Phys. Rev. E} {\bf 49}, 5131 (1993)

\bibitem{vLS93}
M. E. Van Leeuwen and B. Smit, {\em Phys. Rev. Lett.}{\bf 71}, 3991 (1993)

\bibitem{WP92}
D. Wei and G. Patey, {\em Phys. Rev. Lett.}{\bf 68}, 2043 (1992)

\bibitem{WLZ92}
J. J. Weis, D. Levesque and G. J. Zarragoicoechea, {\em Phys. Rev. Lett.}
{\bf 69}, 913 (1992)

\bibitem{WP92a}
D. Wei and G. Patey, {\em Phys. Rev. A}{\bf 46}, 7783 (1992)

\bibitem{WL93}
J. J. Weis and D. Levesque, {\em Phys. Rev. E}{\bf 48}, 3728 (1992)

\bibitem{Joslin93}
C. Joslin and S. Goldman, {\em Mol. Phys.} {\bf 79}, 499 (1993)

\bibitem{Zhang94}
H. Zhang, {\em thesis} (Carnegie Mellon University, 1994)

\bibitem{Huke00}
B. Huke and M. Lucke, {\em Phys. Rev. E}{\bf 62}, 6875 (2000)

\bibitem{WZ92}
M. Widom and H. Zhang, {\em Mat. Res. Soc. Symp.} {\bf 248}, 235 (1992)

\bibitem{Davis99}
S. W. Davis, W. McCausland, H. C. McGahagan, C. T. Tanaka and M. Widom,
{\em Phys. Rev. E}{\bf 59}, 2424 (1999)

\bibitem{Tomanek95}
P. Jund, S. G. Kim, D. Tomanek and J. Hetherington, {\em Phys. Rev. Lett.}
{\bf 74}, 3049 (1995)

\end{thebibliography}

\begin{figure}
\includegraphics[width=5in]{config}
\caption{\label{config} (left) Typical configuration at $\ts=0.2$,
$\rs=0.0125$ projected into $yz$ plane.  Black and white dots indicate
north and south magnetic poles.}
\end{figure}
\newpage
\begin{figure}
\includegraphics[width=5in]{chiz}
\includegraphics[width=5in]{chixy}
\caption{\label{data} Susceptibility data 
versus $\rs$ at fixed $\ts=0.2$.  Top graph illustrates intrinsic
susceptibility $\chi_{int}=\chi_z$, bottom graph illustrates $\chi_{xy}$.
Solid, dotted and dashed lines, respectively, are first, second and
third order Virial series in the strong coupling limit.}
\end{figure}

\end{document}
