Proving spheres are orthogonal 
Given two spheres in $\mathbb{R}^3$:
$x^2+y^2+z^2=2ax; \ \  \ x^2+y^2+z^2 = 2by$
and $a,b>0$, and $\gamma$ the intersection of the spheres, show that for any $p_0 \in \gamma$, the spheres are orthogonal at $p_0$.

I am not sure I fullly understand the problem but here is what I tried:
First of all assuming $p_0$ is some point s.t each sphere can be represented as a function $z_1, z_2 : \mathbb{R}^2 \to \mathbb{R}$ in a neighborhood of $p_0$ (if it isn't we can use another variable), what I think we want to show is that the tangant spaces to $z_1, z_2$ are orthogonal at this point, i.e $\langle\nabla z_1(p_0), \nabla z_2(p_0)\rangle=0.$ But when I calculate this I get that it does not equal 0:
$z_1 = \sqrt{2ax-x^2-y^2}, z_2 = \sqrt{2by-x^2-y^2}$
$\nabla z_1 = (\frac{2a-2x}{2 \sqrt{2ax-x^2-y^2}}, \frac{-2y}{2 \sqrt{2ax-x^2-y^2}})$
$\nabla z_2 = (\frac{-2x}{2 \sqrt{2by-x^2-y^2}}, \frac{2b-2y}{2 \sqrt{2by-x^2-y^2}})$
$\langle\nabla z_1(p_0), \nabla z_2(p_0)\rangle \neq 0.$
Can someone explain where is my mistake?
 A: Let$$f(x,y,z)=(x-a)^2+y^2+z^2\quad\text{and}\quad g(x,y,z)=x^2+(y-b)^2+z^2.$$Then your spheres are the surfaces$$\{(x,y,z)\in\Bbb R^3\mid f(x,y,z)=a^2\}\quad\text{and}\quad\{(x,y,z)\in\Bbb R^3\mid g(x,y,z)=b^2\}.$$If $p=(x_0,y_0,z_0)\in\gamma$ then the tangent plane to the first sphere at $p$ is orthogonal to$$\nabla f(x_0,y_0,z_0)=\bigl(2(x_0-a),2y_0,2z_0\bigr)$$and the tangent plane to the second sphere at $p$ is orthogonal to$$\nabla g(x_0,y_0,z_0)=\bigl(2x_0,2(y_0-b),2z_0\bigr).$$And we have\begin{align}\bigl\langle\nabla f(x_0,y_0,z_0),\nabla g(x_0,y_0,z_0)\bigr\rangle&=\bigl\langle\bigl(2(x_0-a),2y_0,2z_0\bigr),\bigl(2x_0,2(y_0-b),2z_0\bigr)\bigr\rangle\\&=4x_0^{\,2}+4y_0^{\,2}+4z_0^{\,2}-2ax_0-2by_0\\&=2(x_0^{\,2}+y_0^{\,2}+z_0^{\,2}-ax_0)+2(x_0^{\,2}+y_0^{\,2}+z_0^{\,2}-by_0)\\&=0+0\text{ (since $p\in\gamma$)}\\&=0.\end{align}And, since the vectors $\nabla f(x_0,y_0,z_0)$ and $\nabla g(x_0,y_0,z_0)$ are orthogonal, then so are tangent planes at $p$, since one of them is orthogonal to the first vector, whereas the other one is orthogonal to the second vector.
A: The gradients are not applied properly.
Let $f(x,y,z)=x^2+y^2+z^2-2ax$ and $g(x,y,z)=x^2+y^2+z^2-2by$.
Then the spheres are given by $f(x,y,z)=0$ respectively $g(x,y,z)=0$.
The gradients of $f$ and $g$ are orthogonal to their tangent spaces. So to verify if the tangent spaces are orthogonal, we need that for every $(x,y,z)$ on both spheres that $\langle \nabla f(x,y,z),\nabla g(x,y,z)\rangle=0$.
A: First solution:
Here is a quick analytic solution first.
Let $p_0=(x_0,y_0,z_0)$ be a point on both spheres.
Then the tangent plane in $p_0$ to the one or the other sphere has the equation:
$$
\begin{aligned}
2x_0(x-x_0) + 2y_0(y-y_0) + 2z_0(z-z_0) - 2a(x-x_0) &=0\ ,\\
2x_0(x-x_0) + 2y_0(y-y_0) + 2z_0(z-z_0) - 2b(y-y_0) &=0\ ,
\end{aligned}
$$
so that two (unnormalized) normal vectors to the two planes are $(x_0-a,y_0,z_0)$ and $(x_0, y_0-b,z_0)$. Their scalar product is:
$$
x_0(x_0-a)+y_0(y_0-b)+z_0^2=\frac 12\Big( \ (x_0^2+y_0^2+z_0^2-2ax_0)+(x_0^2+y_0^2+z_0^2-2by_0)\ \Big)=0\ .
$$

Second solution:
Trying to go as the OP goes is also possible, even if complicated. The locally defined function $z_1$ is only a part of the parametrization
$$
(x,y)\to(x,y,z_1(x,y))\ .
$$
Letting $x$ variate around $x_0$, and keeping $y=y_0$ we obtain a tangent vector $v_x$ to the sphere, similarly keeping $x=x_0$ and letting $y$ variate around $y_0$ we get an other tangent vector $v_y$. These two vectors, and the normal vector $n_1=v_x\times v_y$,  are:
$$
\begin{aligned}
v_x &= (1,0,\nabla_x z_1(x_0,y_0)) =\left(1, 0, \frac{a-x_0}{\sqrt{2ax_0-x_0^2-y_0^2}}\right)\ ,\\
v_y &= (0,1,\nabla_y z_1(x_0,y_0)) =\left(0, 1, \frac{0-x_0}{\sqrt{2ax_0-x_0^2-y_0^2}}\right)\ ,\\
n_1 &= v_x\times v_y 
=
\left(
\ -\frac{a-x_0}{\sqrt{2ax_0-x_0^2-y_0^2}}\ ,
\ -\frac{0-y_0}{\sqrt{2ax_0-x_0^2-y_0^2}}\ ,
\ 1\
\right)
\\
&\sim \left(\ x_0-a\ ,\ y_0\ ,\ \sqrt{2ax_0-x_0^2-y_0^2}\ \right)
\\
&=(\ x_0-a\ ,\ y_0\ ,\ z_0\ )\ .
\end{aligned}
$$
The computation of the corresponding normal vector using the second parametrization delivers analogously the value
$$
(\ x_0\ ,\ y_0-b\ ,\ z_0\ )\ .
$$
As in the first solution, these vectors are orthogonal.

Third solution:
A geometrical solution is as follows.
The first sphere is a sphere $(S_a)$ centered in $(a,0,0)$ with radius $a$.
The second sphere is a sphere $S_b$ centered in $(0,b,0)$ with radius $b$.
The inversion with center $O=(0,0,0)$ and power $2ab$ transforms $S_a$ in the plane with equation $x=b$, and $S_b$ in the plane with equation $y=a$. Obviously, these planes intersect orthogonally, the inversion is an orthogonal transformation, we apply it again, so $S_a$ and $S_b$ intersect orthogonally.
