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 [itex]\Phi[/itex] does not depend on the azimuthal angle [itex]\phi[/itex] and is also an even function of [itex]z[/itex]. It also satisfies the Laplace equation [itex]\nabla^{2} \, \Phi = 0[/itex] everywhere outside the spheres with the boundary conditions that [itex]\Phi \rightarrow 0[/itex] at infinity and [itex]\Phi = \Phi_{0}[/itex] 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 [itex]a[/itex], then, the equation for such a sphere in cylindrical coordinates is:
[tex]
\rho^{2} + (z - a)^{2} = a^{2}[/tex]
[tex]
\rho^{2} + z^{2} - 2 a z = 0[/tex]
The differential equation for this family of surfaces is obtained by eliminating the parameter [itex]a[/itex]:
[tex]
2 \rho + 2 z z' - 2 a z' = 0 \Rightarrow 2 z (\rho + z z') - 2 a z z' = 0[/tex]
[tex]
z' (\rho^{2} + z^{2}) - 2 z (\rho + z z') = 0[/tex]
[tex]
(\rho^{2} - z^{2}) z' - 2 \rho z = 0[/tex]
A family of orthogonal surfaces is obtaned by the substitution [itex]z \rightarrow -1/z[/itex]:
[tex]
(\rho^{2} - z^{2}) (-\frac{1}{z'}) - 2 \rho z = 0[/tex]
[tex]
\rho^{2} - z^{2} + 2 \rho z z' = 0[/tex]
[tex]
z' = \frac{z^{2} - \rho^{2}}{2 \rho z}[/tex]
This is a homogeneous differential equation and the variables can be separated by applying the substitution [itex]u = z/\rho \Rightarrow z = \rho u, z' = u + \rho u'[/itex]:
[tex]
\rho u' + u = \frac{u^{2} - 1}{2 u}[/tex]
[tex]
\rho u' = - \frac{u^{2} + 1}{2 u}[/tex]
[tex]
\frac{2 u du}{u^{2} + 1} + \frac{d\rho}{\rho} = 0[/tex]
[tex]
\ln{(u^{2} + 1)} + \ln{\rho} = c'[/tex]
[tex]
\rho (u^{2} + 1) = C[/tex]
[tex]
\rho^{2} + z^{2} - C \rho = 0[/tex]
But, this determines a sphere of radius [itex]|b|, \; b = C/2[/itex] lying on the [itex]\rho[/itex]-axis and touching the origin.
Let us use a and b as new coordinates instead of [itex]\rho[/itex] and z. We need to solve the system of equations:
[tex]
\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}}[/tex]
With these definitions, the metric becomes:
[tex]
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}[/tex]
and the Laplace operator takes the form:
[tex]
\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[/tex]
Assuming azimuthal symmetry, this potential does not depend on [itex]\phi[/itex]. However, the equation does not allow separation of the variables a and b.