  We consider the Casimir interaction between two spheres in
$(D+1)$-dimensional Minkowski spacetime due to the vacuum fluctuations of
scalar fields. We consider combinations of Dirichlet and Neumann boundary
conditions. The TGTG formula of the Casimir interaction energy is derived. The
computations of the T matrices of the two spheres are straightforward. To
compute the two G matrices, known as translation matrices, which relate the
hyper-spherical waves in two spherical coordinate frames differ by a
translation, we generalize the operator approach employed in [IEEE Trans.
Antennas Propag. \textbf{36}, 1078 (1988)]. The result is expressed in terms of
an integral over Gegenbauer polynomials. Using our expression for the Casimir
interaction energy, we derive the large separation and small separation
asymptotic expansions of the Casimir interaction energy. In the large
separation regime, we find that the Casimir interaction energy is of order
$L^{-2D+3}$, $L^{-2D+1}$ and $L^{-2D-1}$ respectively for Dirichlet-Dirichlet,
Dirichlet-Neumann and Neumann-Neumann boundary conditions, where $L$ is the
center-to-center distance of the two spheres. In the small separation regime,
we confirm that the leading term of the Casimir interaction agrees with the
proximity force approximation, which is of order $d^{-\frac{D+1}{2}}$, where
$d$ is the distance between the two spheres. Another main result of this work
is the analytic computations of the next-to-leading order term in the small
separation asymptotic expansion. This term is computed using careful order
analysis as well as perturbation method. We find that when $D$ is large, the
ratio of the next-to-leading order term to the leading order term is linear in
$D$, indicating a larger correction at higher dimensions.
