Proof of orthogonality of unit vectors in orthogonal curvilinear coordinate system Context
Suppose $\mathbf{R}(q_1,q_2,q_3)=\mathbf{r}(x,y,z)$ represents a position vector in physical space in a curvilinear coordinate system defined by $q_i=q_i(x,y,z)$ for $i=1,2,3$. The reverse mapping $x=x(q_1,q_2,q_3)$ (similarly for $y$ and $z$) also exists.
The curvilinear system is orthogonal if, at the intersection point of the planes $q_i=\text{constant}$, the normals of these planes are mutually perpendicular. So, at any point $(x_0,y_0,z_0)$
\begin{align}
\nabla q_1\cdot\nabla q_2=0=\nabla q_2\cdot\nabla q_3=\nabla q_3\cdot\nabla q_1. \tag{1}
\end{align}
The unit vectors ($\mathbf{e}_i$) are defined as
\begin{align}
\mathbf{u}_i=&\ \frac{\partial \mathbf{r}}{\partial q_i},\\
\mathbf{e}_i=\frac{\mathbf{u}_i}{\lVert\mathbf{u}_i\rVert}=&\ \frac{\partial \mathbf{r}/\partial q_i}{\lVert\partial \mathbf{r}/\partial q_i\rVert}.
\end{align}
Question
Prove that $\mathbf{e}_i\cdot\mathbf{e}_j=\delta_{ij}$ (Kronecker Delta).
Note
One can either prove

*

*$\mathbf{u}_i\times\nabla q_i=\mathbf{0}$ (i.e.; $\mathbf{u}_i$ is parallel to $\nabla q_i$) so that orthogonality follows from eq. (1),

*Or, prove that $\mathbf{u}_i\cdot\nabla q_j=0$ for $i\neq j$
Using method 2, one way could be to consider
$$\text{d}q_1=1\,\text{d}q_1+0\,\text{d}q_2+0\,\text{d}q_3=\frac{\partial q_1}{\partial x}\text{d}x + \frac{\partial q_1}{\partial y}\text{d}y + \frac{\partial q_1}{\partial z}\text{d}z. \tag{2}$$
Expanding $\text{d}x$ as
$$\text{d}x=\frac{\partial x}{\partial q_1}\text{d}q_1 + \frac{\partial x}{\partial q_2}\text{d}q_2 + \frac{\partial x}{\partial q_3}\text{d}q_3,$$
and substituting back in eq. (2) along with similar expansions for $\text{d}y$ and $\text{d}z$, it can be shown by equating coefficients of $\text{d}q_i$ in eq. (2) that
$$1=\frac{\partial q_1}{\partial x}\frac{\partial x}{\partial q_1} + \frac{\partial q_1}{\partial y}\frac{\partial y}{\partial q_1} + \frac{\partial q_1}{\partial z}\frac{\partial z}{\partial q_1}=\nabla q_1\cdot\mathbf{u}_1,\\
0=\frac{\partial q_1}{\partial x}\frac{\partial x}{\partial q_2} + \frac{\partial q_1}{\partial y}\frac{\partial y}{\partial q_2} + \frac{\partial q_1}{\partial z}\frac{\partial z}{\partial q_2}=\nabla q_1\cdot\mathbf{u}_2,\\
0=\frac{\partial q_1}{\partial x}\frac{\partial x}{\partial q_3} + \frac{\partial q_1}{\partial y}\frac{\partial y}{\partial q_3} + \frac{\partial q_1}{\partial z}\frac{\partial z}{\partial q_3}=\nabla q_1\cdot\mathbf{u}_3.$$
Similarly proceeding with $\text{d}q_2$ and $\text{d}q_3$ proves the result. But this analysis could be done for any curvilinear system. Where was assumption of orthogonal curvilinear system (eq. (1)) invoked? What is wrong in this proof?
Edit
Thanks to Kenny Wong's answer for pointing out that

*

*The relation $\mathbf{u}_i\cdot\nabla q_j=\delta_{ij}$ holds true in any curvilinear coordinate system (need not be orthogonal). There is nothing wrong with the proof shown in the question. It only proves $\mathbf{u}_i\cdot\nabla q_j=\delta_{ij}$ but not that $\mathbf{u}_i\cdot\mathbf{u}_j=0$ for $i\neq j$.

*For an orthogonal curvilinear coordinate system, this result ($\mathbf{u}_i\cdot\nabla q_j=\delta_{ij}$) along with eq. (1) proves that $\mathbf{u}_i$ is parallel to $\nabla q_i$ and hence $\mathbf{u}_i\cdot\mathbf{u}_j=0$ for $i\neq j$.

 A: The source of the confusion is that $\mathbf{ u}_i $ and $\nabla q_i $ are not the same thing. They are not necessarily even parallel!
Let me give you an example. We'll work in two dimensions, with
$$ q_1 = x + y, \ \ \ q_2 = y.$$
Inverting this, we get
$$ x = q_1 - q_2, \ \ \ y = q_2.$$
So
$$ \mathbf{u}_1 = (1, 0), \ \ \ \mathbf{u}_2 = (-1, 1),$$
$$ \nabla q_1 = (1, 1), \ \ \ \nabla q_2 = (0, 1).$$
As you can see the $\mathbf{u}_i$'s are not the same as the $\nabla q_i$'s, not even up to scaling!
It is true that
$$ \mathbf{u}_i . \nabla q_j = \delta_{ij},$$
as you correctly proved in your question.
But you're interested in something different: you're interested in whether or not $\mathbf{u}_i . \mathbf u_j = 0$ for $i \neq j$. And since the $\mathbf{u}_i$'s are not parallel to the $\nabla q_i$'s, we simply can't tell whether or not $\mathbf{u}_i . \mathbf u_j = 0$ for $i \neq j$, knowing only that $\mathbf{u}_i . \nabla q_j = \delta_{ij}$.
I believe people call $\{ \mathbf{u}_1, \mathbf{u}_2 \}$ the "original" basis and $\{ \nabla q_1, \nabla q_2 \}$ the "dual" basis. They are two different bases, as this example shows.

Returning to the original question, I think the thing to do is to show that $\mathbf{u}_i$ and $\nabla q_i$ are the same, up to scaling, on the condition that the coordinates are orthogonal.
This is actually quite easy to do. Let's write
$$ \mathbf{u}_1 = c_1 \nabla q_1 + c_2 \nabla q_2,$$
where $c_1, c_2$ are cooefficients that we must determine.
If we dot both sides with $\nabla q_1$, or with $\nabla q_2$, we get
$$ 1 = \mathbf{u}_1 . \nabla q_1 = c_1 \nabla q_1 . \nabla q_1 + c_2 \nabla q_2 . \nabla q_1  = c_1 \| \nabla q_1\|^2$$
$$ 0 =  \mathbf{u}_1 . \nabla q_2 = c_1 \nabla q_1 . \nabla q_2 + c_2 \nabla q_2 . \nabla q_2  = c_2 \| \nabla q_2\|^2 $$
Let me clarify what I've just done:


*

*On the left hand side, I used the relations $\mathbf{u_i} . \nabla q_j = \delta_{ij}$ which you proved in generality for all coordinate systems.

*On the right hand side, I used the fact that $\nabla q_1 . \nabla q_2 = 0$, which is true only for an orthogonal choice of coodinates.
Anyway, this exercise tells us that $c_1 = 1 / \| \nabla q_1 \|^2$ and $c_2 = 0$, so
$$ \mathbf{u}_1 = \frac{\nabla q_1}{\| \nabla q_1 \|^2}.$$
And if we do the same for $\mathbf{u}_2$, we get
$$ \mathbf{u}_2 = \frac{\nabla q_2}{\| \nabla q_2 \|^2}.$$
So, under the assumption that the coordinates are orthogonal, we find that the two sets of basis vectors  $\{ \mathbf{u}_1, \mathbf{u}_2 \}$ and $\{\nabla q_1, \nabla q_2 \}$ are the same thing, up to rescaling. Since $\{\nabla q_1, \nabla q_2 \}$  is an orthogonal basis by assumption, it must follow that $\{ \mathbf{u}_1, \mathbf{u}_2 \}$ is orthogonal too.
