This is a very hard problem.
EDIT:
Let's take the z-axis along the line connecting the two centers of the two spheres, and the point where the spheres touch each other lie in the xy-plane (it follows that this point coincides with the origin).
From the geometry of the problem, it follows then that the z-axis is an axis of symmetry and that the xy-plane is a plane of symmetry for the problem. Therefore, the potential \Phi does not depend on the azimuthal angle \phi and is also an even function of z. It also satisfies the Laplace equation \nabla^{2} \, \Phi = 0 everywhere outside the spheres with the boundary conditions that \Phi \rightarrow 0 at infinity and \Phi = \Phi_{0} on the spheres.
Let us try and solve this boundary value problem. Before we go any further, let us try and "invent" a new coordinate system, where one set of coordinate surfaces are spheres with the center on the z-axis and such that the origin always lies on the sphere. If the sphere has radius a, then, the equation for such a sphere in cylindrical coordinates is:
<br />
\rho^{2} + (z - a)^{2} = a^{2}<br />
<br />
\rho^{2} + z^{2} - 2 a z = 0<br />
The differential equation for this family of surfaces is obtained by eliminating the parameter a:
<br />
2 \rho + 2 z z' - 2 a z' = 0 \Rightarrow 2 z (\rho + z z') - 2 a z z' = 0<br />
<br />
z' (\rho^{2} + z^{2}) - 2 z (\rho + z z') = 0<br />
<br />
(\rho^{2} - z^{2}) z' - 2 \rho z = 0<br />
A family of orthogonal surfaces is obtaned by the substitution z \rightarrow -1/z:
<br />
(\rho^{2} - z^{2}) (-\frac{1}{z'}) - 2 \rho z = 0<br />
<br />
\rho^{2} - z^{2} + 2 \rho z z' = 0<br />
<br />
z' = \frac{z^{2} - \rho^{2}}{2 \rho z}<br />
This is a homogeneous differential equation and the variables can be separated by applying the substitution u = z/\rho \Rightarrow z = \rho u, z' = u + \rho u':
<br />
\rho u' + u = \frac{u^{2} - 1}{2 u}<br />
<br />
\rho u' = - \frac{u^{2} + 1}{2 u}<br />
<br />
\frac{2 u du}{u^{2} + 1} + \frac{d\rho}{\rho} = 0<br />
<br />
\ln{(u^{2} + 1)} + \ln{\rho} = c'<br />
<br />
\rho (u^{2} + 1) = C<br />
<br />
\rho^{2} + z^{2} - C \rho = 0<br />
But, this determines a sphere of radius |b|, \; b = C/2 lying on the \rho-axis and touching the origin.
Let us use a and b as new coordinates instead of \rho and z. We need to solve the system of equations:
<br />
\left\{\begin{array}{lcl}<br />
\rho^{2} + z^{2} - 2 a z & = & 0 \\<br />
<br />
\rho^{2} + z^{2} - 2 b \rho & = & 0<br />
\end{array}\right. \Rightarrow \left\{\begin{array}{l}<br />
\rho = \frac{2 a^{2} b}{a^{2} + b^{2}} \\<br />
<br />
z = \frac{2 a b^{2}}{a^{2} + b^{2}}<br />
\end{array}\right. \Rightarrow r^{2} = \rho^{2} + z^{2} = \frac{4 a^{2} b^{2}}{a^{2} + b^{2}}<br />
With these definitions, the metric becomes:
<br />
ds^{2} = d\rho^{2} + dz^{2} + \rho^{2} d\phi^{2} = \frac{4(b^{4} \, da^{2} + a^{4} \, db^{2}}{(a^{2} + b^{2})^{2}} + \frac{4 a^{4} \, b^{2}}{(a^{2} + b^{2})^{2}} d\phi^{2}<br />
and the Laplace operator takes the form:
<br />
\nabla^{2} \Phi = \frac{<br />
(a^{2} + b^{2})^{3}}{4 b^{4}} \, \frac{\partial}{\partial a}\left(\frac{a^{4}}{a^{2} + b^{2}} \, \frac{\partial \Phi}{\partial a}\right) + \frac{(a^{2} + b^{2})^{3}}{4 a^{2} \, b^{3}} \, \frac{\partial}{\partial b}\left(\frac{b}{a^{2} + b^{2}} \, \frac{\partial \Phi}{\partial b}\right) + \frac{(a^{2} + b^{2})^{2}}{4 a^{4} b^{2}} \, \frac{\partial^{2} \Phi}{\partial \phi^{2}} = 0<br />
Assuming azimuthal symmetry, this potential does not depend on \phi. However, the equation does not allow separation of the variables a and b.
)? What does that make V in the capacitance formula C=\frac{Q}{V}?
