  We have developed an exact, general method to compute Casimir interactions
between a finite number of compact objects of arbitrary shape and separation.
Here, we present details of the method for a scalar field to illustrate our
approach in its most simple form; the generalization to electromagnetic fields
is outlined in Ref. [1]. The interaction between the objects is attributed to
quantum fluctuations of source distributions on their surfaces, which we
decompose in terms of multipoles. A functional integral over the effective
action of multipoles gives the resulting interaction. Each object's shape and
boundary conditions enter the effective action only through its scattering
matrix. Their relative positions enter through universal translation matrices
that depend only on field type and spatial dimension. The distinction of our
method from the pairwise summation of two-body potentials is elucidated in
terms of the scattering processes between three objects. To illustrate the
power of the technique, we consider Robin boundary conditions $\phi -\lambda
\partial_n \phi=0$, which interpolate between Dirichlet and Neumann cases as
$\lambda$ is varied. We obtain the interaction between two such spheres
analytically in a large separation expansion, and numerically for all
separations. The cases of unequal radii and unequal $\lambda$ are studied. We
find sign changes in the force as a function of separation in certain ranges of
$\lambda$ and see deviations from the proximity force approximation even at
short separations, most notably for Neumann boundary conditions.
