Divergence in curvilinear coordinates

In a curvilinear system with coordinates: $u_1, u_1, u_3$.

I want to understand why: $$\nabla\cdot \vec{V} \neq \frac{1}{h_1} \frac{\partial V_1}{\partial u_1} +\frac{1}{h_2} \frac{\partial V_2}{\partial u_2} + \frac{1}{h_3}\frac{\partial V_3}{\partial u_3}. \tag1$$

My reasoning to explain why eq.(1) is not true is the following:

Taking the first term of the dot product:

$$\frac{\hat{e_1} }{h_1}\frac{\partial}{\partial u_1} (V_1 \hat{e_1}) = \frac{\hat{e_1}}{h_1} \big( \hat{e_1} \frac{\partial V_1}{\partial u_1} + V_1 \frac{\partial \hat{e_1}}{\partial u_1} \big)$$

This seems rare to me since the definition of the dot product is:

$$A \cdot B = A_1B_1 + A_2B_2 + A_3B_3 $$

And using that definition the divergence should be eq.(1) (with an equality).

I would appreciate any references that explain this, since all I can find is the formula for the divergence in curvilinear systems but without a deduction.

Defining a "del" vector $\nabla=\mathbf{e}_i\frac{\partial}{\partial x^i}$ (summation on repeated indices assumed) is sort-of a mistake. This isn't a vector (field). It is just a nice notational device, because in this case, in cartesian coordinates, we have $$ \text{grad}(f)=\mathbf{e}_i\frac{\partial f}{\partial x^ i},$$ which looks like $\nabla f$, where $\nabla$ is the "vector" defined above and we have $$ \text{div}(\mathbf{A})=\frac{\partial A^ i}{\partial x^i}, $$ which sort of looks like $\nabla\cdot\mathbf{A}$, etc. But this is just a useful way to remember these formulas.

Now let's note several things. We can define a matrix whose $i,j$-th element is $\partial A^ i/\partial x^ j$. Let's call this matrix $\nabla\otimes\mathbf{A}$. We clearly have $$ \text{div}(\mathbf{A})=\text{Tr}(\nabla\otimes\mathbf{A})=\sum_i\frac{\partial A^ i}{\partial x^ i}. $$ (I'm gonna suspend the automatic summation from now on.)

If we have an orthogonal curvilinear coordinate system $(u^1,u^2,u^3)$, we can define the following quantities:

• Coordinate basis vectors $$ \mathbf{g}_i=\frac{\partial \mathbf{r}}{\partial u^ i}, $$

• "Reciprocal" coordinate basis vectors $$ \mathbf{g}^i=\nabla u^ i, $$

• An orthonormal frame $$ \hat{\mathbf{e}}_i=\frac{1}{h_i}\mathbf{g}_i, $$ where $h_i=\sqrt{\mathbf{g}_i\cdot\mathbf{g}_i}$.

These quantities have the property that $$ \mathbf{g}^ i\cdot\mathbf{g}_j=\delta^ i_j, $$ since ($x^ i$ are cartesian coordinates and $\mathbf{e}_i$ are cartesian basis vectors) we have $$ \mathbf{g}^i\cdot\mathbf{g}_j=\nabla u^i\cdot\frac{\partial\mathbf{r}}{\partial u^ j}=\sum_k\frac{\partial u^ i}{\partial x^k}\mathbf{e}_k\cdot\sum_l\frac{\partial x^l}{\partial u^ j}\mathbf{e}_l=\sum_{kl}\frac{\partial u^ i}{\partial x^k}\frac{\partial x^l}{\partial u_j}\delta_{kl}=\sum_k \frac{\partial u^ i}{\partial x^k}\frac{\partial x^k}{\partial u_j}=\delta^i_j.$$

We further introduce a matrix $g_{ij}$ given as $g_{ij}=\mathbf{g}_i\cdot\mathbf{g}_j$. Because of orthogonality, we have $g_{ij}=h_i^2\delta_{ij}$. We also introduce $g^{ij}=\mathbf{g}^i\cdot\mathbf{g}^j$, feel free to check that $g^{ij}=\frac{1}{h_i^2}\delta_{ij}$, thus $g_{ij}$ and $g^{ij}$ are inverse matrices.

If $\mathbf{A}$ is a vector field, we associate three kinds of components with $\mathbf{A}$. We can have

• $\mathbf{A}=\sum_i A^ i\mathbf{g}_i$, and the $A^ i$ are called contravariant components,

• $\mathbf{A}=\sum_i A_i\mathbf{g}^i$, and the $A_i$ are called covariant components and

• $\mathbf{A}=\sum_i \hat{A}_i\hat{\mathbf{e}}_i$ and the $\hat{A}_i$ are called physical components (here whether you put the index upwards or downwards doesn't matter!).

It is clear that we have $A^ i=\mathbf{g}^i\cdot\mathbf{A}$ and $A_i=\mathbf{g}_i\cdot\mathbf{A}$, and what is less clear is that differentiation with respect to the coordinates $\partial/\partial u^ i$ produce covariant components, eg. we have for any scalar function $f$, $\nabla f=\sum_i\mathbf{g}^i\frac{\partial f}{\partial u^ i}$ (check the definition of $\mathbf{g}^i$!).

We define Christoffel symbols of the second kind as $$ \Gamma^ k_{ij}=\frac{1}{2}\sum_l g^{kl}\left(\frac{\partial g_{jl}}{\partial u^ i}+\frac{\partial g_{il}}{\partial u^ j}-\frac{\partial g_{ij}}{\partial u^ l}\right). $$ Check that $\Gamma^k_{ij}=\mathbf{g}^k\cdot\frac{\partial \mathbf{g}_i}{\partial u^ j}$ (I really don't want to derive this identity right now).

In curvilinear coordinates, the basis vectors also depend on positions, so every time you differentiate a vector field, you need to make sure to take the variation of the basis vectors also into account, so we calculate the divergence as $$ \text{div}(\mathbf{A})=\sum_i\mathbf{g}^i\cdot\frac{\partial\mathbf{A}}{\partial u^i}=\sum_{ij}\mathbf{g}^i\left(\frac{\partial A^j}{\partial u^ i}\mathbf{g}_j+A^j\frac{\partial \mathbf{g}_j}{\partial u^ i}\right)=\sum_{ij}\left(\frac{\partial A^ i}{\partial u^i}+\Gamma^i_{ji}A^j\right). $$

We now calculate $$\Gamma^i_{ji}=\sum\frac{1}{2}g^{il}\left(\frac{\partial g_{jl}}{\partial u^i}+\frac{\partial g_{il}}{\partial u^j}-\frac{\partial g_{ij}}{\partial u^l}\right)=\sum\frac{1}{2}g^{il}\frac{\partial g_{il}}{\partial u^j},$$ which, in matrix notation ($(g_{ij})=g$) can be written as $$ \Gamma^i_{ji}=\frac{1}{2}\text{Tr}\left(g^{-1}\frac{\partial g}{\partial u^j}\right), $$ which, using Jacobi's formula, can be written as $$ \sum_i\Gamma^ i_{ji}=\frac{1}{2}\frac{1}{\det g}\frac{\partial}{\partial u^ j}\det g=\frac{1}{2}\frac{\partial}{\partial u^ j}\ln\det g=\frac{\partial}{\partial u^j}\ln\sqrt{\det g}=\frac{1}{\sqrt{\det g}}\frac{\partial}{\partial u^j}\sqrt{\det g}. $$

Plugging this back into the divergence formula gives $$\text{div}(\mathbf{A})=\sum_{ij}\left(\frac{\partial A^ i}{\partial u^i}+\Gamma^i_{ji}A^j\right)=\sum_{ij}\left(\frac{\partial A^i}{\partial u^i}+\frac{1}{\sqrt{\det g}}\frac{\partial}{\partial u^j}\sqrt{\det g}A^j\right)= \\ =\sum_{i}\frac{1}{\sqrt{\det g}}\frac{\partial}{\partial u^i}(\sqrt{\det g}A^i). $$

But since $g=\text{diag}(h_1^2,h_2^2,h_3^2)$

, we have $\det g=(h_1h_2h_3)^2$, so $$ \text{div}(\mathbf{A})=\frac{1}{h_1h_2h_3}\left\{\frac{\partial}{\partial u^1}(h_1h_2h_3A^1)+\frac{\partial}{\partial u^2}(h_1h_2h_3A^2)+\frac{\partial}{\partial u^3}(h_1h_2h_3A^3)\right\}= \\ =\frac{1}{h_1h_2h_3}\left\{\frac{\partial}{\partial u^1}(h_2h_3\hat{A}_1)+\frac{\partial}{\partial u^2}(h_1h_3\hat{A}_2)+\frac{\partial}{\partial u^3}(h_1h_2\hat{A}_3)\right\}, $$

which is the formula for divergence in terms of the physical components.