It seems that the potential inside your surface is indeed not constant.

I faintly remember that a charge distrubition on a surface with varying curvature is not constant, but I am not sure.

Look at ( I exchanged x and z )

xx = {Cos[phi] Sin[th], Sin[phi] Sin[th], 2 Cos[th]};

xx[[1]]^2 + xx[[2]]^2 + xx[[3]]^2/4 // Simplify

ParametricPlot3D[xx, {phi, 0, 2 Pi}, {th, 0, Pi}]

and

e1 = D[xx, th];
e2 = D[xx, phi];
df1 = Cross[e1, e2];
df2 = FullSimplify[Sqrt[df1.df1]]

dd1 = xx - {x, y, z};
dd2 = Sqrt[dd1.dd1] // FullSimplify

Then

pot[x_, y_, z_] :=
 NIntegrate[Sqrt[(5 - 3 Cos[2 th]) Sin[th]^2]/Sqrt[2]/
  Sqrt[(z - 2 Cos[th])^2 + (x - Cos[phi] Sin[th])^2 + (y - 
     Sin[phi] Sin[th])^2],
  {phi, 0, 2 Pi}, {th, 0, Pi}]

Obviously the potential calculated in this manner is not constant:

ppot = Table[pot[x, 0, 0], {x, -2.5, 2.5, .09}];
ListLinePlot[ppot]

Sorry. The above method is simply wrong.

I think you cannot calculate a surface integral in a manner like this. I think this is due to that NIntegrate choses a lot of points in the region of integration (discretize the region) and almost none fulfills the condition in the Boole-statement. So the result is rubbish.

Try to calculate the surface of a (unit) sphere

NIntegrate[ Boole[x^2 + y^2 + z^2 == 1.0], {x, -2, 2}, {y, -2, 2}, {z, -2, 2}]

The result is wrong.

But when calculating the volume a sufficent amount of points fulfilling the Boole-condition is found and you get a correct result

NIntegrate[Boole[x^2 + y^2 + z^2 <= 1.0], {x, -2, 2}, {y, -2, 2}, {z, -2, 2}]
N[4 Pi /3]

So forget my proposal.