How to Calculate Divergence in Curvilinear Coordinates

Divergence measures the net outward flux of a vector field per unit volume. In Cartesian coordinates, its familiar form is

$$\nabla\cdot\mathbf A=\frac{\partial A_x}{\partial x}+\frac{\partial A_y}{\partial y}+\frac{\partial A_z}{\partial z}.$$

In curvilinear coordinates, coordinate surfaces may have different spacings and different areas. The scale factors that describe this geometry must therefore appear in the formula. This article gives the general formula, the especially useful orthogonal-coordinate version, and worked cylindrical and spherical examples.

The general coordinate formula

Let $(q^1,q^2,q^3)$ be coordinates, and let $g_{ij}$ be the metric tensor. Write

$$g=\det(g_{ij}).$$

If the vector field is expressed using contravariant components $A^i$, then its divergence is

$$\boxed{\nabla\cdot\mathbf A =\frac{1}{\sqrt{g}}\frac{\partial}{\partial q^i}\left(\sqrt{g}\,A^i\right)}$$

where repeated indices are summed over $i=1,2,3$. Written without index summation notation, this is

$$\nabla\cdot\mathbf A =\frac{1}{\sqrt{g}} \left[ \frac{\partial}{\partial q^1}(\sqrt{g}A^1) +\frac{\partial}{\partial q^2}(\sqrt{g}A^2) +\frac{\partial}{\partial q^3}(\sqrt{g}A^3) \right].$$

The components $A^i$ here are contravariant coordinate components. They are not always the same as the components measured along unit vectors. This distinction is the source of many errors in curvilinear-coordinate calculations.

Orthogonal coordinates and scale factors

Most physics and vector-calculus applications use orthogonal coordinates. Let the line element be

$$ds^2=h_1^2(dq^1)^2+h_2^2(dq^2)^2+h_3^2(dq^3)^2,$$

where $h_1,h_2,h_3$ are the scale factors. In this case,

$$\sqrt{g}=h_1h_2h_3.$$

If $A_1,A_2,A_3$ denote the physical components of the vector field along the corresponding local unit vectors, then

$$\boxed{ \nabla\cdot\mathbf A =\frac{1}{h_1h_2h_3} \left[ \frac{\partial}{\partial q^1}(h_2h_3A_1) +\frac{\partial}{\partial q^2}(h_3h_1A_2) +\frac{\partial}{\partial q^3}(h_1h_2A_3) \right] }$$

This is the standard divergence formula for an orthogonal curvilinear coordinate system.

The physical and contravariant components are related by

$$A^i=\frac{A_i}{h_i}.$$

Thus, if you use the scale-factor formula above, use physical components $A_i$ consistently.

Cylindrical coordinates

Cylindrical coordinates $(r,\phi,z)$ are defined by

$$x=r\cos\phi,\qquad y=r\sin\phi,\qquad z=z.$$

Their scale factors are

$$h_r=1,\qquad h_\phi=r,\qquad h_z=1.$$

For

$$\mathbf A=A_r\,\mathbf e_r+A_\phi\,\mathbf e_\phi+A_z\,\mathbf e_z,$$

the divergence is therefore

$$\boxed{ \nabla\cdot\mathbf A =\frac{1}{r}\frac{\partial}{\partial r}(rA_r) +\frac{1}{r}\frac{\partial A_\phi}{\partial\phi} +\frac{\partial A_z}{\partial z} }$$

Cylindrical example

Take

$$\mathbf A=r^2\,\mathbf e_r+r\sin\phi\,\mathbf e_\phi+z^2\,\mathbf e_z.$$

Thus $A_r=r^2$, $A_\phi=r\sin\phi$, and $A_z=z^2$. Substitute into the formula:

$$\begin{aligned} \nabla\cdot\mathbf A &=\frac{1}{r}\frac{\partial}{\partial r}(r\cdot r^2) +\frac{1}{r}\frac{\partial}{\partial\phi}(r\sin\phi) +\frac{\partial}{\partial z}(z^2)\\ &=\frac{1}{r}\frac{\partial}{\partial r}(r^3) +\frac{1}{r}(r\cos\phi)+2z\\ &=3r+\cos\phi+2z. \end{aligned}$$

The answer is

$$\boxed{\nabla\cdot\mathbf A=3r+\cos\phi+2z.}$$

Spherical coordinates

Use the convention $(r,\theta,\phi)$ in which $r$ is radius, $\theta$ is the polar angle measured from the positive $z$-axis, and $\phi$ is the azimuthal angle in the $xy$-plane:

$$x=r\sin\theta\cos\phi,\qquad y=r\sin\theta\sin\phi,\qquad z=r\cos\theta.$$

The scale factors are

$$h_r=1,\qquad h_\theta=r,\qquad h_\phi=r\sin\theta.$$

For

$$\mathbf A=A_r\,\mathbf e_r+A_\theta\,\mathbf e_\theta+A_\phi\,\mathbf e_\phi,$$

the divergence becomes

$$\boxed{ \nabla\cdot\mathbf A =\frac{1}{r^2}\frac{\partial}{\partial r}(r^2A_r) +\frac{1}{r\sin\theta}\frac{\partial}{\partial\theta}(\sin\theta A_\theta) +\frac{1}{r\sin\theta}\frac{\partial A_\phi}{\partial\phi} }$$

Spherical example

Consider the radial field

$$\mathbf A=r^2\,\mathbf e_r.$$

Here $A_r=r^2$ and $A_\theta=A_\phi=0$. Therefore

$$\begin{aligned} \nabla\cdot\mathbf A &=\frac{1}{r^2}\frac{\partial}{\partial r}(r^2r^2)\\ &=\frac{1}{r^2}\frac{\partial}{\partial r}(r^4)\\ &=\frac{1}{r^2}(4r^3)=4r. \end{aligned}$$

Hence

$$\boxed{\nabla\cdot(r^2\mathbf e_r)=4r.}$$

More generally, for a purely radial spherical field $\mathbf A=f(r)\mathbf e_r$,

$$\nabla\cdot\mathbf A =\frac{1}{r^2}\frac{d}{dr}\left(r^2f(r)\right).$$

For example, $\mathbf A=\mathbf e_r/r^2$ has zero divergence wherever $r\ne0$, although it has a singular source at the origin in the distributional sense.

Common mistakes

Differentiating components without scale factors

The expression $\partial A_r/\partial r+\partial A_\theta/\partial\theta+\partial A_\phi/\partial\phi$ is not the divergence in spherical coordinates. The coordinate directions and area elements change with position.

Confusing angle conventions

Some books call the azimuthal angle $\theta$ and the polar angle $\phi$, reversing the names used above. The formula must follow the definitions of the coordinates, not merely their symbols.

Mixing component types

The orthogonal formula uses components along unit vectors, such as $A_r,A_\theta,A_\phi$. The metric-determinant formula uses contravariant components $A^i$. Substituting one type into the other without the appropriate scale factors gives an incorrect result.

Ignoring singular coordinate locations

Cylindrical formulas contain $1/r$, and spherical formulas contain $1/(r\sin\theta)$. These factors reflect coordinate singularities at the axis or origin. A formula valid away from those locations may need separate analysis there, especially for fields with point or line sources.