Curvilinear coordinates

Differential operators in orthogonal coordinate systems.

In this chapter 2

Tick a unit to add it to your PDF basket.

Orthogonal curvilinear coordinates
- Orthogonal curvilinear coordinates
- Cylindrical and spherical coordinates

Orthogonal curvilinear coordinates

Method

Orthogonal curvilinear coordinates

Adapted

Add to PDF

Builds the orthonormal tangent frame and the scale factors of a curvilinear coordinate system, gives the line, area and volume elements they produce, and derives the four general formulas for grad, div, curl and the Laplacian, with a fully worked example in parabolic cylindrical coordinates.

Prerequisites: Divergence, curl and the Laplacian, Gradient and directional derivative, Vector differential identities, Multiple integrals and change of variables, Implicit and inverse function theorems

Definition(Curvilinear coordinates, scale factors, the unit basis)AdaptedAdapted from Riley, Hobson & Bence (2006), §10.10, pp. 364–365; Arfken, Weber & Harris (2013), §3.10, pp. 183–184, Eqs. (3.128)–(3.132)

Let $D \subseteq \mathbb{R}^3$ be open and let

\mathbf{r} : D \to \mathbb{R}^3, \qquad (u_1, u_2, u_3) \mapsto \mathbf{r}(u_1,u_2,u_3) = \bigl(x(u_1,u_2,u_3),\, y(\cdots),\, z(\cdots)\bigr)

be of class $C^2$ , injective, with a $C^1$ inverse — a curvilinear coordinate system on the image $\mathbf{r}(D)$ . (Riley, Hobson and Bence ask only for continuity, differentiability and a single-valued inverse; $C^2$ is assumed here because the proof below differentiates twice.) Both of our sources allow the conditions to break down on a few exceptional points and curves, such as the origin or the polar axis; here those sets are simply removed from $D$ .

The coordinate surfaces through a point are the three level sets $u_i = \text{const}$ ; each pair of them meets in a coordinate curve, along which only one coordinate varies.
The tangent vectors are $\mathbf{e}_i = \partial\mathbf{r}/\partial u_i$ , tangent to the $u_i$ -curve in the direction of increasing $u_i$ . We assume $\mathbf{e}_i \neq \mathbf{0}$ throughout $D$ .
The scale factors are their lengths, $h_i = \lvert \mathbf{e}_i \rvert = \left\lvert \frac{\partial\mathbf{r}}{\partial u_i} \right\rvert = \left[ \left(\frac{\partial x}{\partial u_i}\right)^{2} + \left(\frac{\partial y}{\partial u_i}\right)^{2} + \left(\frac{\partial z}{\partial u_i}\right)^{2} \right]^{1/2} > 0 ,$ and the unit basis is $\hat{\mathbf{e}}_i = \mathbf{e}_i / h_i$ . A coordinate $u_i$ need not be a length; $h_i$ carries whatever dimension makes $h_i\,du_i$ one.
The system is orthogonal when $\mathbf{e}_i \cdot \mathbf{e}_j = 0$ for $i \neq j$ at every point of $D$ , so that $\hat{\mathbf{e}}_1, \hat{\mathbf{e}}_2, \hat{\mathbf{e}}_3$ is an orthonormal basis at each point. Equivalently, in Arfken's notation the metric coefficients $g_{ij} = \mathbf{e}_i \cdot \mathbf{e}_j$ form a diagonal matrix, $g_{ij} = h_i^2\,\delta_{ij}$ . Throughout this unit the labelling is also chosen so that the frame is right-handed, $\hat{\mathbf{e}}_1 \times \hat{\mathbf{e}}_2 = \hat{\mathbf{e}}_3$ .

A vector field is written in the unit basis as $\mathbf{a} = a_1\hat{\mathbf{e}}_1 + a_2\hat{\mathbf{e}}_2 + a_3\hat{\mathbf{e}}_3$ ; the $a_i$ are its physical components, and they carry the dimensions of $\mathbf{a}$ itself. The frame rotates from point to point, so the $\hat{\mathbf{e}}_i$ must be differentiated along with the $a_i$ whenever an operator is applied to $\mathbf{a}$ .

Method(Line, area and volume elements)AdaptedAdapted from Riley, Hobson & Bence (2006), §10.10, Eqs. (10.55)–(10.58), pp. 365–366; Arfken, Weber & Harris (2013), §3.10, Eqs. (3.133)–(3.137), pp. 184–185

Let $u_1,u_2,u_3$ be an orthogonal curvilinear system on $D$ as above. Then, at every point of $D$ :

Displacement. $\displaystyle d\mathbf{r} = \sum_{i=1}^{3} \frac{\partial\mathbf{r}}{\partial u_i}\,du_i = h_1\,du_1\,\hat{\mathbf{e}}_1 + h_2\,du_2\,\hat{\mathbf{e}}_2 + h_3\,du_3\,\hat{\mathbf{e}}_3 .$
Arc length. Since the $\hat{\mathbf{e}}_i$ are orthonormal, $(ds)^2 = d\mathbf{r}\cdot d\mathbf{r} = h_1^2\,(du_1)^2 + h_2^2\,(du_2)^2 + h_3^2\,(du_3)^2 ,$ with no cross terms; a line integral of $\mathbf{a}$ along a curve is $\int \mathbf{a}\cdot d\mathbf{r} = \sum_i \int a_i\,h_i\,du_i$ .
Area. On a coordinate surface $u_1 = \text{const}$ , the patch spanned by $du_2$ and $du_3$ has area $h_2h_3\,du_2\,du_3$ and unit normal $\pm\hat{\mathbf{e}}_1$ , so the flux of $\mathbf{a}$ through it is $a_1\,h_2h_3\,du_2\,du_3$ ; the two other faces are obtained by cycling $1\to2\to3\to1$ .
Volume. The infinitesimal box spanned by $h_i\,du_i\,\hat{\mathbf{e}}_i$ has volume $dV = \bigl\lvert h_1\,du_1\,\hat{\mathbf{e}}_1 \cdot (h_2\,du_2\,\hat{\mathbf{e}}_2 \times h_3\,du_3\,\hat{\mathbf{e}}_3) \bigr\rvert = h_1h_2h_3\,du_1\,du_2\,du_3 .$ The scalar triple product on the left is the determinant whose columns are the $\mathbf{e}_i$ , so this says that the Jacobian factor of the map $\mathbf{r}$ is $\lvert\det D\mathbf{r}\rvert = h_1h_2h_3$ (change of variables) — orthogonality is what lets the determinant be read off from the three lengths alone.

Method(The four operators in orthogonal curvilinear coordinates)AdaptedAdapted from Riley, Hobson & Bence (2006), §10.10, Eqs. (10.60)–(10.62) and Table 10.4, pp. 367–369; Arfken, Weber & Harris (2013), §3.10, Eqs. (3.138)–(3.143), pp. 185–187

With the hypotheses of the definition, let $\Phi$ be a $C^1$ scalar field (a $C^2$ one for the Laplacian) and $\mathbf{a} = a_1\hat{\mathbf{e}}_1 + a_2\hat{\mathbf{e}}_2 + a_3\hat{\mathbf{e}}_3$ a $C^1$ vector field on $\mathbf{r}(D)$ . Writing $H = h_1h_2h_3$ :

\nabla\Phi = \frac{1}{h_1}\frac{\partial\Phi}{\partial u_1}\,\hat{\mathbf{e}}_1 + \frac{1}{h_2}\frac{\partial\Phi}{\partial u_2}\,\hat{\mathbf{e}}_2 + \frac{1}{h_3}\frac{\partial\Phi}{\partial u_3}\,\hat{\mathbf{e}}_3 ,

\nabla\cdot\mathbf{a} = \frac{1}{H}\left[ \frac{\partial}{\partial u_1}\bigl(h_2h_3\,a_1\bigr) + \frac{\partial}{\partial u_2}\bigl(h_3h_1\,a_2\bigr) + \frac{\partial}{\partial u_3}\bigl(h_1h_2\,a_3\bigr) \right] ,

\nabla\times\mathbf{a} = \frac{1}{H} \begin{vmatrix} h_1\hat{\mathbf{e}}_1 & h_2\hat{\mathbf{e}}_2 & h_3\hat{\mathbf{e}}_3 \\[2pt] \dfrac{\partial}{\partial u_1} & \dfrac{\partial}{\partial u_2} & \dfrac{\partial}{\partial u_3} \\[6pt] h_1a_1 & h_2a_2 & h_3a_3 \end{vmatrix} ,

\nabla^2\Phi = \frac{1}{H}\left[ \frac{\partial}{\partial u_1}\!\left( \frac{h_2h_3}{h_1}\frac{\partial\Phi}{\partial u_1} \right) + \frac{\partial}{\partial u_2}\!\left( \frac{h_3h_1}{h_2}\frac{\partial\Phi}{\partial u_2} \right) + \frac{\partial}{\partial u_3}\!\left( \frac{h_1h_2}{h_3}\frac{\partial\Phi}{\partial u_3} \right) \right] .

The determinant is expanded along its first row, the entries of the middle row acting as operators on the entries of the last; explicitly its $\hat{\mathbf{e}}_1$ component is $\frac{1}{h_2h_3}\bigl[\partial_{u_2}(h_3a_3) - \partial_{u_3}(h_2a_2)\bigr]$ , and the other two follow by cycling. Setting all $h_i = 1$ recovers the Cartesian formulas of Divergence, curl and the Laplacian.

Idea of the proofAdaptedAdapted from Riley, Hobson & Bence (2006), §10.10, pp. 367–369

Only the gradient has to be computed from scratch, and it comes straight from comparing two expressions for the same differential $d\Phi$ : the chain rule gives it in terms of $\partial\Phi/\partial u_i$ , and the definition of the gradient gives it as $\nabla\Phi\cdot d\mathbf{r}$ , where $d\mathbf{r}$ carries the scale factors. Applying that one formula to $\Phi = u_i$ yields the key relation $\nabla u_i = \hat{\mathbf{e}}_i/h_i$ , which turns each basis vector into a gradient: $\hat{\mathbf{e}}_1 = h_1\nabla u_1$ and also $\hat{\mathbf{e}}_1 = \hat{\mathbf{e}}_2\times\hat{\mathbf{e}}_3 = h_2h_3\,\nabla u_2\times\nabla u_3$ . The divergence and the curl then follow from the two vanishing identities $\nabla\cdot(\nabla\phi\times\nabla\psi) = 0$ and $\nabla\times\nabla\phi = \mathbf{0}$ : each kills the awkward term of a product rule and leaves a single gradient, which the first formula evaluates. The Laplacian is the divergence formula fed with the gradient formula.

ProofAdaptedAdapted from Riley, Hobson & Bence (2006), §10.10, pp. 367–369; Arfken, Weber & Harris (2013), §3.10, Eqs. (3.138) and (3.142), pp. 185–186

Throughout, the frame is orthonormal and right-handed, so $\hat{\mathbf{e}}_1 = \hat{\mathbf{e}}_2\times\hat{\mathbf{e}}_3$ , $\hat{\mathbf{e}}_2 = \hat{\mathbf{e}}_3\times\hat{\mathbf{e}}_1$ , $\hat{\mathbf{e}}_3 = \hat{\mathbf{e}}_1\times\hat{\mathbf{e}}_2$ , and the component of any vector along $\hat{\mathbf{e}}_i$ is its dot product with $\hat{\mathbf{e}}_i$ .

1. The gradient. Fix a point and let the coordinates change by $du_1,du_2,du_3$ . The chain rule for a $C^1$ function of three variables gives $d\Phi = \sum_i (\partial\Phi/\partial u_i)\,du_i$ . On the other hand $\Phi$ is a field on $\mathbf{r}(D)$ , so by the defining property of the gradient, $d\Phi = \nabla\Phi\cdot d\mathbf{r}$ . Expand both factors in the unit basis: writing $\nabla\Phi = \sum_i g_i\hat{\mathbf{e}}_i$ and using $d\mathbf{r} = \sum_i h_i\,du_i\,\hat{\mathbf{e}}_i$ from the element statement, orthonormality kills every cross term and leaves $\nabla\Phi\cdot d\mathbf{r} = \sum_i g_i h_i\,du_i$ . The three increments $du_i$ are independent, so the coefficients must agree one by one: $g_i h_i = \partial\Phi/\partial u_i$ . Since $h_i > 0$ we may divide, which is the gradient formula.

2. The reciprocal relation. First, the coordinate functions are twice differentiable as functions of position. The columns of $D\mathbf{r}$ are the three tangent vectors $\mathbf{e}_i$ , which are mutually orthogonal and have lengths $h_i > 0$ , hence linearly independent; so $\lvert\det D\mathbf{r}\rvert = h_1h_2h_3 \neq 0$ at every point of $D$ . Since $\mathbf{r}$ is $C^2$ with nowhere-vanishing Jacobian determinant, the inverse function theorem (Implicit and inverse function theorems) makes it a local $C^2$ diffeomorphism; as it is also injective, $\mathbf{r}(D)$ is open and $\mathbf{r}^{-1}$ is $C^2$ on it. In particular each coordinate function $u_j$ is a $C^2$ function of $(x,y,z)$ — a stronger conclusion than the $C^1$ inverse assumed in the definition, and one we need in steps 3 and 5.

Now apply step 1 to $\Phi = u_j$ . Its partial derivatives are $\partial u_j/\partial u_i = \delta_{ij}$ , so only one term survives:

\nabla u_j = \frac{1}{h_j}\,\hat{\mathbf{e}}_j .

Thus $\hat{\mathbf{e}}_j = h_j \nabla u_j$ : in an orthogonal system the frame tangent to the coordinate curves and the frame normal to the coordinate surfaces coincide.

3. The divergence. Treat $\mathbf{a} = \sum_i a_i\hat{\mathbf{e}}_i$ one term at a time; the divergence is linear. For the first term, use right-handedness and then step 2 twice:

a_1\hat{\mathbf{e}}_1 = a_1\,(\hat{\mathbf{e}}_2\times\hat{\mathbf{e}}_3) = a_1\,(h_2\nabla u_2)\times(h_3\nabla u_3) = \bigl(a_1h_2h_3\bigr)\,\bigl(\nabla u_2\times\nabla u_3\bigr).

Apply the product rule $\nabla\cdot(f\mathbf{B}) = \nabla f\cdot\mathbf{B} + f\,\nabla\cdot\mathbf{B}$ with $f = a_1h_2h_3$ and $\mathbf{B} = \nabla u_2\times\nabla u_3$ . The second term vanishes, because $\nabla\cdot(\nabla\phi\times\nabla\psi) = 0$ for any two $C^2$ scalars (Vector differential identities); this is where $\mathbf{r}$ being $C^2$ is used, since by step 2 it is what makes $u_2$ and $u_3$ twice differentiable functions of position. Hence, using step 2 once more to rewrite $\mathbf{B} = \hat{\mathbf{e}}_2\times\hat{\mathbf{e}}_3/(h_2h_3) = \hat{\mathbf{e}}_1/(h_2h_3)$ ,

\nabla\cdot(a_1\hat{\mathbf{e}}_1) = \nabla\bigl(a_1h_2h_3\bigr)\cdot\frac{\hat{\mathbf{e}}_1}{h_2h_3} .

The dot product with $\hat{\mathbf{e}}_1$ picks out the first component of the gradient, which step 1 gives as $(1/h_1)\partial_{u_1}(a_1h_2h_3)$ . So $\nabla\cdot(a_1\hat{\mathbf{e}}_1) = H^{-1}\partial_{u_1}(h_2h_3a_1)$ . Cycling the labels $1\to2\to3\to1$ — which preserves right-handedness — gives the other two terms, and adding them gives the divergence formula.

4. The Laplacian. By definition $\nabla^2\Phi = \nabla\cdot(\nabla\Phi)$ . Step 1 says that the physical components of $\nabla\Phi$ are $a_i = (1/h_i)\,\partial\Phi/\partial u_i$ , which are $C^1$ when $\Phi$ is $C^2$ . Substituting them into step 3 gives the stated formula at once. Note that no mixed derivative $\partial^2/\partial u_i\partial u_j$ with $i \neq j$ can appear: orthogonality is exactly what removes them.

5. The curl. Again work one term at a time, now using $\hat{\mathbf{e}}_1 = h_1\nabla u_1$ from step 2, so that $a_1\hat{\mathbf{e}}_1 = (a_1h_1)\nabla u_1$ . The product rule $\nabla\times(f\mathbf{B}) = \nabla f\times\mathbf{B} + f\,\nabla\times\mathbf{B}$ with $f = a_1h_1$ and $\mathbf{B} = \nabla u_1$ loses its second term because $\nabla\times\nabla\phi = \mathbf{0}$ for a $C^2$ scalar (Vector differential identities), leaving

\nabla\times(a_1\hat{\mathbf{e}}_1) = \nabla\bigl(a_1h_1\bigr)\times\frac{\hat{\mathbf{e}}_1}{h_1} .

Expand the gradient by step 1. Its $\hat{\mathbf{e}}_1$ part crosses with $\hat{\mathbf{e}}_1$ to zero, and the other two use $\hat{\mathbf{e}}_2\times\hat{\mathbf{e}}_1 = -\hat{\mathbf{e}}_3$ and $\hat{\mathbf{e}}_3\times\hat{\mathbf{e}}_1 = \hat{\mathbf{e}}_2$ :

\nabla\times(a_1\hat{\mathbf{e}}_1) = \frac{1}{h_3h_1}\frac{\partial (h_1a_1)}{\partial u_3}\,\hat{\mathbf{e}}_2 \;-\; \frac{1}{h_1h_2}\frac{\partial (h_1a_1)}{\partial u_2}\,\hat{\mathbf{e}}_3 .

Cycling gives the contributions of $a_2\hat{\mathbf{e}}_2$ and $a_3\hat{\mathbf{e}}_3$ . Collect the three results: the terms along $\hat{\mathbf{e}}_1$ come from $a_3$ and $a_2$ and add up to

\frac{1}{h_2h_3}\left[ \frac{\partial (h_3a_3)}{\partial u_2} - \frac{\partial (h_2a_2)}{\partial u_3} \right] ,

and cyclically for the other two. Multiplying and dividing by $h_1$ makes this $H^{-1}h_1\bigl[\partial_{u_2}(h_3a_3) - \partial_{u_3}(h_2a_2)\bigr]$ , which is precisely the cofactor of $h_1\hat{\mathbf{e}}_1$ in the stated determinant; the same holds for the other two entries, so the determinant is a compact way of writing the three components at once. $\square$

Remark(Where the hypotheses are used)Original

Orthogonality enters three times: in step 1, to make $\nabla\Phi\cdot d\mathbf{r}$ separate into three independent terms; in step 3, through $\hat{\mathbf{e}}_1 = \hat{\mathbf{e}}_2\times\hat{\mathbf{e}}_3$ , which needs the frame to be orthonormal and right-handed, not merely orthogonal; and in the volume element, to make the determinant equal $h_1h_2h_3$ . Drop orthogonality and none of the four formulas survives: a non-orthogonal system has $g_{ij} \neq 0$ off the diagonal, the Laplacian acquires mixed derivatives $g^{ij}\partial_i\partial_j$ , and the right language is tensor analysis with raised and lowered indices.
$h_i > 0$ , that is $\partial\mathbf{r}/\partial u_i \neq \mathbf{0}$ , is needed to divide by $h_i$ , to define $\hat{\mathbf{e}}_i$ at all, and (in step 2) to keep $\det D\mathbf{r}$ away from zero. It is exactly the condition that fails on the excluded sets: at $r = 0$ in spherical polars $h_\theta = r$ vanishes, and on the polar axis $h_\varphi = r\sin\theta$ vanishes. Equivalently $H = h_1h_2h_3$ is the Jacobian factor, and the coordinates degenerate where it vanishes.
Right-handed labelling is a convenience in the proof but a real hypothesis only for the curl. Step 3 writes $\hat{\mathbf{e}}_1 = \hat{\mathbf{e}}_2\times\hat{\mathbf{e}}_3$ , which is false in a left-handed frame; but there the minus sign appears twice — once in $a_1\hat{\mathbf{e}}_1 = -(a_1h_2h_3)\,\nabla u_2\times\nabla u_3$ and once in $\mathbf{B} = \nabla u_2\times\nabla u_3 = -\hat{\mathbf{e}}_1/(h_2h_3)$ — and the two cancel in $\nabla f\cdot\mathbf{B}$ . So the gradient, divergence and Laplacian formulas hold whatever the handedness. The curl does not: swapping two coordinates, say $u_1 \leftrightarrow u_2$ , reverses the sign of the determinant and hence of $\nabla\times\mathbf{a}$ . Spherical polars $(r,\theta,\varphi)$ and cylindrical polars $(\rho,\varphi,z)$ are both listed in right-handed order.
The frame is not constant. Because $\hat{\mathbf{e}}_i$ varies, $\nabla\cdot\hat{\mathbf{e}}_1$ and $\nabla\times\hat{\mathbf{e}}_1$ need not vanish even though $\hat{\mathbf{e}}_1$ has constant length. For the same reason the vector Laplacian is not obtained by applying $\nabla^2$ to each physical component separately; see Cylindrical and spherical coordinates.

Remark(How to read the formulas)AdaptedAdapted from Riley, Hobson & Bence (2006), §10.10, closing paragraph, p. 369; Arfken, Weber & Harris (2013), §3.10, Eq. (3.140), p. 186, and Exercise 3.10.4(a), p. 197

The non-orthogonal case. Riley, Hobson and Bence close §10.10 by pointing the reader to the general tensor treatment, which is where a system with $g_{ij} \neq 0$ off the diagonal has to be handled.
What the divergence formula means. The combination $h_2h_3\,a_1$ is the flux of $\mathbf{a}$ per unit $du_2\,du_3$ through a face $u_1 = \text{const}$ , and $H\,du_1du_2du_3$ is the volume of the cell. So the formula is the net outflow per unit volume, computed face by face, and that is how Arfken derives it: the two opposite faces of a cell contribute $\pm h_2h_3a_1$ evaluated at $u_1 \pm \tfrac12 du_1$ , and it is essential that $h_2h_3$ as well as $a_1$ be evaluated on each face, since the two faces have different areas. The curl formula reads the same way with circulations around the faces instead of fluxes through them.
A worked consequence. Putting $a_1 = 1$ , $a_2 = a_3 = 0$ in the divergence formula gives $\nabla\cdot\hat{\mathbf{e}}_1 = H^{-1}\partial_{u_1}(h_2h_3)$ , which is Arfken's Exercise 3.10.4(a) and quantifies the previous remark's point that the frame is not constant.

NoteAdaptedAdapted from Lee (2013), Problems 16-11 and 16-13, p. 436

The same formulas in Riemannian language. For a Riemannian metric $g$ in arbitrary smooth local coordinates $x^i$ , Lee's Problem 16-11 states the divergence of a vector field with contravariant components $X^i$ as

\operatorname{div} X = \frac{1}{\sqrt{\det g}}\,\frac{\partial}{\partial x^i}\Bigl( X^i \sqrt{\det g} \Bigr) ,

and his Problem 16-13 states the Laplacian as $\Delta u = -\frac{1}{\sqrt{\det g}}\partial_i\bigl( g^{ij}\sqrt{\det g}\,\partial_j u \bigr)$ . The minus sign is Lee's convention: he defines the geometric Laplacian as $\Delta u = -\operatorname{div}(\operatorname{grad} u)$ , the negative of the operator $\nabla^2$ used here and throughout this handbook, so that its eigenvalues are non-negative; he warns explicitly that conventions differ between authors. Up to that sign the two statements contain ours. For an orthogonal system, $g_{ij} = h_i^2\delta_{ij}$ gives $\sqrt{\det g} = h_1h_2h_3 = H$ and $g^{ij} = h_i^{-2}\delta_{ij}$ , while a physical component $a_i$ is related to the contravariant one by $X^i = a_i/h_i$ (the basis $\mathbf{e}_i$ has length $h_i$ ). Substituting turns Lee's two formulas into the divergence and Laplacian stated above.

Method(Setting up a new orthogonal coordinate system)Original

Given a candidate map $\mathbf{r}(u_1,u_2,u_3)$ :

Differentiate. Compute the three tangent vectors $\mathbf{e}_i = \partial\mathbf{r}/\partial u_i$ in Cartesian components. Everything else is algebra on these.
Test orthogonality by checking $\mathbf{e}_1\cdot\mathbf{e}_2 = \mathbf{e}_2\cdot\mathbf{e}_3 = \mathbf{e}_3\cdot\mathbf{e}_1 = 0$ identically. If any product is non-zero the formulas of this unit do not apply.
Read off the scale factors $h_i = \lvert\mathbf{e}_i\rvert$ and record where they vanish; delete those sets from $D$ and state the surviving domain, together with the ranges of the $u_i$ on which $\mathbf{r}$ is injective.
Fix the orientation. Check $\hat{\mathbf{e}}_1\times\hat{\mathbf{e}}_2 = +\hat{\mathbf{e}}_3$ ; if it comes out $-\hat{\mathbf{e}}_3$ , swap two of the coordinates before using the curl.
Write down the elements: $(ds)^2$ , the three face areas $h_2h_3\,du_2du_3$ and its cyclic partners, and $dV = H\,du_1du_2du_3$ . These are what integrals need.
Substitute into the four formulas. Simplify $H$ , $h_2h_3/h_1$ and its cyclic partners first; in most useful systems these collapse, which is what makes the Laplacian tractable.
Check against a known field. Take a simple Cartesian field, express its physical components in the new frame, apply the formulas, and compare with the Cartesian answer. Two or three such tests catch almost every algebraic slip.

How it fails. Step 2 is the one that is skipped and should not be: elliptic, parabolic and bipolar systems are orthogonal, but many natural-looking maps (a sheared or a body-fitted grid) are not, and applying these formulas to them is simply wrong. Step 3 is the second: a formula that is correct on $D$ says nothing on the degenerate set, and a field that looks singular on the axis in these coordinates may be perfectly smooth there.

Example(Parabolic cylindrical coordinates)Original

Define, for $u > 0$ and $v > 0$ and all real $z$ ,

x = \tfrac12\bigl(u^2 - v^2\bigr), \qquad y = uv, \qquad z = z .

Geometry and injectivity. In complex notation $x + iy = \tfrac12 (u + iv)^2$ , so the map is the squaring map $w \mapsto w^2/2$ applied to the open first quadrant of the $w$ -plane, whose image is the open upper half plane $y > 0$ ; squaring is injective there because it doubles arguments lying in $(0,\pi/2)$ . The surfaces $u = \text{const}$ and $v = \text{const}$ are confocal parabolic cylinders with the $z$ -axis as common focal line, opening in the $-x$ and $+x$ directions.
Tangent vectors. $\mathbf{e}_u = (u, v, 0)$ , $\mathbf{e}_v = (-v, u, 0)$ , $\mathbf{e}_z = (0,0,1)$ .
Orthogonality. $\mathbf{e}_u\cdot\mathbf{e}_v = -uv + vu = 0$ , and both are perpendicular to $\mathbf{e}_z$ since their third components vanish. The system is orthogonal.
Scale factors. $h_u = h_v = \sqrt{u^2+v^2}$ and $h_z = 1$ , so $H = u^2 + v^2 > 0$ on the whole domain. The two transverse scale factors are equal — the map is conformal in the $xy$ -plane, as the complex form predicts, with $\lvert dw^2/2 \,/\, dw\rvert = \lvert w\rvert = \sqrt{u^2+v^2}$ .
Orientation. $\hat{\mathbf{e}}_u \times \hat{\mathbf{e}}_v = (u,v,0)\times(-v,u,0)/(u^2+v^2) = (0,0,u^2+v^2)/(u^2+v^2) = \hat{\mathbf{e}}_z$ , so $(u,v,z)$ is already right-handed.
Elements. $(ds)^2 = (u^2+v^2)\bigl[(du)^2 + (dv)^2\bigr] + (dz)^2$ and $dV = (u^2+v^2)\,du\,dv\,dz$ .
The Laplacian. Here $h_vh_z/h_u = 1$ , $h_zh_u/h_v = 1$ and $h_uh_v/h_z = u^2+v^2$ , so the general formula collapses to $\nabla^2\Phi = \frac{1}{u^2+v^2}\left( \frac{\partial^2\Phi}{\partial u^2} + \frac{\partial^2\Phi}{\partial v^2} \right) + \frac{\partial^2\Phi}{\partial z^2} .$
First check (step 7 of the method). Take $\Phi = x^2 + y^2$ , for which $\nabla^2\Phi = 4$ in Cartesian coordinates. In the new variables $x^2 + y^2 = \tfrac14(u^2-v^2)^2 + u^2v^2 = \tfrac14(u^2+v^2)^2$ . Then $\Phi_{uu} = 3u^2 + v^2$ and $\Phi_{vv} = u^2 + 3v^2$ , whose sum is $4(u^2+v^2)$ ; dividing by $u^2+v^2$ gives $\nabla^2\Phi = 4$ , as required.
Second check: a vector field. Take the rigid-rotation field $\mathbf{F} = -y\,\mathbf{i} + x\,\mathbf{j}$ , for which $\nabla\cdot\mathbf{F} = 0$ and $\nabla\times\mathbf{F} = 2\mathbf{k}$ . Its physical components are $F_u = \mathbf{F}\cdot\hat{\mathbf{e}}_u = \frac{-uv\cdot u + \tfrac12(u^2-v^2)\cdot v}{\sqrt{u^2+v^2}} = -\tfrac12 v\sqrt{u^2+v^2}, \qquad F_v = \tfrac12 u\sqrt{u^2+v^2}, \qquad F_z = 0 ,$ so that $h_uF_u = -\tfrac12 v(u^2+v^2)$ and $h_vF_v = \tfrac12 u(u^2+v^2)$ . The divergence formula gives $\nabla\cdot\mathbf{F} = \frac{1}{u^2+v^2}\Bigl[ \partial_u\bigl(-\tfrac12 v(u^2+v^2)\bigr) + \partial_v\bigl(\tfrac12 u(u^2+v^2)\bigr) \Bigr] = \frac{-uv + uv}{u^2+v^2} = 0 .$ For the curl, the $\hat{\mathbf{e}}_u$ and $\hat{\mathbf{e}}_v$ components vanish because $F_z = 0$ and nothing depends on $z$ , while the $\hat{\mathbf{e}}_z$ component is $\frac{1}{h_uh_v}\Bigl[ \partial_u(h_vF_v) - \partial_v(h_uF_u) \Bigr] = \frac{1}{u^2+v^2}\left[ \frac{3u^2+v^2}{2} + \frac{u^2+3v^2}{2} \right] = \frac{2(u^2+v^2)}{u^2+v^2} = 2 .$ So $\nabla\times\mathbf{F} = 2\hat{\mathbf{e}}_z = 2\mathbf{k}$ , matching the Cartesian value. Both checks are independent of the coordinate algebra that produced them, which is the point of step 7.

ApplicationOriginal

Separation of variables. The reason for cataloguing these systems is that the Laplacian above separates in a handful of them, and each separation names a family of special functions. The example's parabolic cylinder coordinates are the case behind diffraction at a straight edge. The edge lies along the $z$ -axis, the problem is two-dimensional, so put $\Phi = U(u)V(v)$ with no $z$ -dependence in the Helmholtz equation $\nabla^2\Phi + k^2\Phi = 0$ . The collapsed Laplacian of the example gives $(u^2+v^2)^{-1}(U''V + UV'') + k^2UV = 0$ ; multiply by $(u^2+v^2)/(UV)$ and the equation becomes

\frac{U''}{U} + k^2u^2 \;+\; \frac{V''}{V} + k^2v^2 \;=\; 0 .

What makes this separate is that $H = u^2+v^2$ is a sum of a function of $u$ and a function of $v$ : once it has been cleared from the denominator, the $u$ 's and the $v$ 's sit in different terms. Setting the first group equal to a constant $\lambda$ and the second to $-\lambda$ ,

U'' + \bigl(k^2u^2 - \lambda\bigr)U = 0, \qquad V'' + \bigl(k^2v^2 + \lambda\bigr)V = 0 ,

two parabolic cylinder equations: the substitution $s = \sqrt{2k}\,e^{i\pi/4}u$ sends the first to Weber's standard form $d^2U/ds^2 + (\nu + \tfrac12 - \tfrac14 s^2)U = 0$ with $\nu + \tfrac12 = i\lambda/(2k)$ .

The straight-edge field is the elementary member of this family, $\lambda = ik$ . Then $U = e^{iku^2/2}$ solves the first equation (indeed $U'' = (ik - k^2u^2)U$ ), and writing $V = e^{-ikv^2/2}f(v)$ reduces the second to $f'' = 2ikv f'$ , so that $f(v) = \int_0^v e^{ik\tau^2}\,d\tau$ — a Fresnel integral. Since $u^2 - v^2 = 2x$ , the product is

\Phi = e^{ikx}\int_0^{v} e^{ik\tau^2}\,d\tau ,

a plane wave modulated by a Fresnel integral of the single parabolic coordinate $v$ , and it vanishes on $v = 0$ , which is exactly the half-plane $y = 0$ , $x > 0$ . The modulus of that integral, traced across the geometric shadow boundary, is the fringe pattern at a straight edge. (If instead a propagating factor $e^{ik_zz}$ is kept, the same algebra goes through with $k^2$ replaced by the transverse wavenumber $k_\perp^2 = k^2 - k_z^2$ ; the separation is only useful when $k_\perp \neq 0$ , since at $k_\perp = 0$ the quadratic terms disappear and $U$ and $V$ are mere exponentials.) Cylindrical and spherical polars give Bessel functions and spherical harmonics in the same way (Cylindrical and spherical coordinates); prolate spheroidal coordinates give the functions used for a dipole antenna, and toroidal coordinates those used for the field of a current ring.

Conservative discretisation. In computational fluid dynamics and plasma simulation the grid is usually curvilinear, following the geometry of a nozzle, a blade passage or a tokamak flux surface. Finite-volume schemes discretise the divergence in exactly the form derived here: the flux through a face is $h_2h_3\,a_1$ times the face's coordinate extent, and the cell volume is $H\,\Delta u_1\Delta u_2\Delta u_3$ . Summing over neighbouring cells then telescopes exactly, because each interior face is counted once with each sign, so the discrete scheme conserves mass and momentum to machine precision — a property that is lost if the operator is expanded out into partial derivatives of the components first.

Reading a field on the axis. The remark's warning about the degenerate set is a practical one. Spacecraft and tokamak diagnostics report vector data in the natural curvilinear frame, where the components of a perfectly smooth field blow up or become undefined where a scale factor vanishes. Interpolating such data near the axis without first transforming back to a frame that is continuous there is a standard source of spurious signals.

Reference table

Cylindrical and spherical coordinates

Adapted

Add to PDF

Tabulates grad, div, curl and the Laplacian, together with the line, area and volume elements, in cylindrical and spherical polar coordinates, and derives the Laplacian twice — from the general curvilinear formula and directly by the chain rule — with worked examples in pipe flow and electrostatics.

Prerequisites: Orthogonal curvilinear coordinates, Divergence, curl and the Laplacian

Definition(Cylindrical polar coordinates)AdaptedAdapted from Riley, Hobson & Bence (2006), §10.9.1, Eqs. (10.44)–(10.52), pp. 357–360; Arfken, Weber & Harris (2013), §3.10, Eqs. (3.144)–(3.146), pp. 188–189

Riley, Hobson and Bence give the ranges $\rho \ge 0$ , $0 \le \varphi < 2\pi$ , $-\infty < z < \infty$ ; here $\rho = 0$ is excluded, for the reason given below. So, for $\rho > 0$ , $0 \le \varphi < 2\pi$ and $-\infty < z < \infty$ ,

x = \rho\cos\varphi, \qquad y = \rho\sin\varphi, \qquad z = z .

The coordinate surfaces are the cylinders $\rho = \text{const}$ about the $z$ -axis, the half-planes $\varphi = \text{const}$ through that axis, and the planes $z = \text{const}$ . Differentiating $\mathbf{r} = \rho\cos\varphi\,\mathbf{i} + \rho\sin\varphi\,\mathbf{j} + z\,\mathbf{k}$ and normalising gives the unit basis

\hat{\mathbf{e}}_\rho = \cos\varphi\,\mathbf{i} + \sin\varphi\,\mathbf{j}, \qquad \hat{\mathbf{e}}_\varphi = -\sin\varphi\,\mathbf{i} + \cos\varphi\,\mathbf{j}, \qquad \hat{\mathbf{e}}_z = \mathbf{k} ,

which is orthonormal and right-handed, with $\hat{\mathbf{e}}_\rho \times \hat{\mathbf{e}}_\varphi = \hat{\mathbf{e}}_z$ . The first two rotate with $\varphi$ and so must be differentiated along with the components of a field. The scale factors are

h_\rho = 1, \qquad h_\varphi = \rho, \qquad h_z = 1, \qquad H = h_\rho h_\varphi h_z = \rho .

The $z$ -axis, where $h_\varphi = 0$ and $\varphi$ is undefined, is excluded. Dropping the $z$ coordinate gives plane polar coordinates $(\rho,\varphi)$ , for which $dA = \rho\,d\rho\,d\varphi$ .

Definition(Spherical polar coordinates)AdaptedAdapted from Riley, Hobson & Bence (2006), §10.9.2, Eqs. (10.53)–(10.54), pp. 361–362; Arfken, Weber & Harris (2013), §3.10, Eqs. (3.152)–(3.153), pp. 190–191

In the convention used here and by both sources, $\theta$ is the polar angle measured from the positive $z$ -axis and $\varphi$ the azimuth measured from the positive $x$ -axis. Arfken gives the ranges $0 \le r < \infty$ , $0 \le \theta \le \pi$ , $0 \le \varphi < 2\pi$ (Riley, Hobson and Bence write $r \ge 0$ for the first); here $r = 0$ is excluded, for the reason given below. So, for $r > 0$ , $0 \le \theta \le \pi$ and $0 \le \varphi < 2\pi$ ,

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

The coordinate surfaces are spheres $r = \text{const}$ centred at the origin, right circular cones $\theta = \text{const}$ with vertex at the origin and axis the $z$ -axis, and half-planes $\varphi = \text{const}$ through that axis. The unit basis is

\begin{aligned} \hat{\mathbf{e}}_r &= \sin\theta\cos\varphi\,\mathbf{i} + \sin\theta\sin\varphi\,\mathbf{j} + \cos\theta\,\mathbf{k}, \\ \hat{\mathbf{e}}_\theta &= \cos\theta\cos\varphi\,\mathbf{i} + \cos\theta\sin\varphi\,\mathbf{j} - \sin\theta\,\mathbf{k}, \\ \hat{\mathbf{e}}_\varphi &= -\sin\varphi\,\mathbf{i} + \cos\varphi\,\mathbf{j} , \end{aligned}

again orthonormal and right-handed, with $\hat{\mathbf{e}}_r \times \hat{\mathbf{e}}_\theta = \hat{\mathbf{e}}_\varphi$ ; all three depend on position. The scale factors are

h_r = 1, \qquad h_\theta = r, \qquad h_\varphi = r\sin\theta, \qquad H = r^2\sin\theta .

Arfken records that $\theta$ and $\varphi$ are both undefined at $r = 0$ and that $\varphi$ is undefined at $\theta = 0$ and $\theta = \pi$ ; these are exactly the places where a scale factor vanishes, and they are excluded.

Reference table(Operators in cylindrical polar coordinates)AdaptedAdapted from Riley, Hobson & Bence (2006), Table 10.2, p. 360; Arfken, Weber & Harris (2013), Eqs. (3.146)–(3.151), pp. 189–190

For a scalar field $\Phi(\rho,\varphi,z)$ and a vector field $\mathbf{a} = a_\rho\hat{\mathbf{e}}_\rho + a_\varphi\hat{\mathbf{e}}_\varphi + a_z\hat{\mathbf{e}}_z$ , off the $z$ -axis:

d\mathbf{r} = d\rho\,\hat{\mathbf{e}}_\rho + \rho\,d\varphi\,\hat{\mathbf{e}}_\varphi + dz\,\hat{\mathbf{e}}_z, \qquad (ds)^2 = (d\rho)^2 + \rho^2(d\varphi)^2 + (dz)^2, \qquad dV = \rho\,d\rho\,d\varphi\,dz ,

and the three coordinate faces carry $\rho\,d\varphi\,dz$ (normal $\hat{\mathbf{e}}_\rho$ ), $d\rho\,dz$ (normal $\hat{\mathbf{e}}_\varphi$ ) and $\rho\,d\rho\,d\varphi$ (normal $\hat{\mathbf{e}}_z$ ).

\nabla\Phi = \frac{\partial\Phi}{\partial\rho}\,\hat{\mathbf{e}}_\rho + \frac{1}{\rho}\frac{\partial\Phi}{\partial\varphi}\,\hat{\mathbf{e}}_\varphi + \frac{\partial\Phi}{\partial z}\,\hat{\mathbf{e}}_z ,

\nabla\cdot\mathbf{a} = \frac{1}{\rho}\frac{\partial}{\partial\rho}\bigl(\rho\,a_\rho\bigr) + \frac{1}{\rho}\frac{\partial a_\varphi}{\partial\varphi} + \frac{\partial a_z}{\partial z} ,

\nabla\times\mathbf{a} = \frac{1}{\rho} \begin{vmatrix} \hat{\mathbf{e}}_\rho & \rho\,\hat{\mathbf{e}}_\varphi & \hat{\mathbf{e}}_z \\[2pt] \dfrac{\partial}{\partial\rho} & \dfrac{\partial}{\partial\varphi} & \dfrac{\partial}{\partial z} \\[6pt] a_\rho & \rho\,a_\varphi & a_z \end{vmatrix} , \qquad \nabla^2\Phi = \frac{1}{\rho}\frac{\partial}{\partial\rho}\!\left( \rho\frac{\partial\Phi}{\partial\rho} \right) + \frac{1}{\rho^2}\frac{\partial^2\Phi}{\partial\varphi^2} + \frac{\partial^2\Phi}{\partial z^2} .

The vector Laplacian is not componentwise. Arfken's Eq. (3.151) gives

(\nabla^2\mathbf{a})_\rho = \nabla^2 a_\rho - \frac{a_\rho}{\rho^2} - \frac{2}{\rho^2}\frac{\partial a_\varphi}{\partial\varphi}, \qquad (\nabla^2\mathbf{a})_\varphi = \nabla^2 a_\varphi - \frac{a_\varphi}{\rho^2} + \frac{2}{\rho^2}\frac{\partial a_\rho}{\partial\varphi}, \qquad (\nabla^2\mathbf{a})_z = \nabla^2 a_z ,

where $\nabla^2$ on the right is the scalar operator above. Only the $z$ -component is the naive one.

Reference table(Operators in spherical polar coordinates)AdaptedAdapted from Riley, Hobson & Bence (2006), Table 10.3, p. 363; Arfken, Weber & Harris (2013), Eqs. (3.154)–(3.160), pp. 191–192

For a scalar field $\Phi(r,\theta,\varphi)$ and a vector field $\mathbf{a} = a_r\hat{\mathbf{e}}_r + a_\theta\hat{\mathbf{e}}_\theta + a_\varphi\hat{\mathbf{e}}_\varphi$ , away from the origin and the polar axis:

d\mathbf{r} = dr\,\hat{\mathbf{e}}_r + r\,d\theta\,\hat{\mathbf{e}}_\theta + r\sin\theta\,d\varphi\,\hat{\mathbf{e}}_\varphi, \qquad (ds)^2 = (dr)^2 + r^2(d\theta)^2 + r^2\sin^2\theta\,(d\varphi)^2 ,

dV = r^2\sin\theta\,dr\,d\theta\,d\varphi, \qquad d\Omega = \sin\theta\,d\theta\,d\varphi \ \text{ with } \int d\Omega = 4\pi ,

the sphere $r = \text{const}$ carrying the area element $r^2\,d\Omega$ with normal $\hat{\mathbf{e}}_r$ .

\nabla\Phi = \frac{\partial\Phi}{\partial r}\,\hat{\mathbf{e}}_r + \frac{1}{r}\frac{\partial\Phi}{\partial\theta}\,\hat{\mathbf{e}}_\theta + \frac{1}{r\sin\theta}\frac{\partial\Phi}{\partial\varphi}\,\hat{\mathbf{e}}_\varphi ,

\nabla\cdot\mathbf{a} = \frac{1}{r^2}\frac{\partial}{\partial r}\bigl(r^2 a_r\bigr) + \frac{1}{r\sin\theta}\frac{\partial}{\partial\theta}\bigl(\sin\theta\,a_\theta\bigr) + \frac{1}{r\sin\theta}\frac{\partial a_\varphi}{\partial\varphi} ,

\nabla\times\mathbf{a} = \frac{1}{r^2\sin\theta} \begin{vmatrix} \hat{\mathbf{e}}_r & r\,\hat{\mathbf{e}}_\theta & r\sin\theta\,\hat{\mathbf{e}}_\varphi \\[2pt] \dfrac{\partial}{\partial r} & \dfrac{\partial}{\partial\theta} & \dfrac{\partial}{\partial\varphi} \\[6pt] a_r & r\,a_\theta & r\sin\theta\,a_\varphi \end{vmatrix} ,

\nabla^2\Phi = \frac{1}{r^2}\frac{\partial}{\partial r}\!\left( r^2\frac{\partial\Phi}{\partial r} \right) + \frac{1}{r^2\sin\theta}\frac{\partial}{\partial\theta}\!\left( \sin\theta\,\frac{\partial\Phi}{\partial\theta} \right) + \frac{1}{r^2\sin^2\theta}\frac{\partial^2\Phi}{\partial\varphi^2} .

Riley, Hobson and Bence note the equivalent radial form, useful for shortening calculations:

\frac{1}{r^2}\frac{\partial}{\partial r}\!\left( r^2\frac{\partial\Phi}{\partial r} \right) = \frac{1}{r}\frac{\partial^2}{\partial r^2}\bigl( r\Phi \bigr) .

Here too the vector Laplacian is not componentwise; Arfken's Eq. (3.160) gives all three components, each of which couples $a_r$ , $a_\theta$ and $a_\varphi$ .

Idea of the proofOriginal

There are two honest ways to get these tables, and it is worth seeing both. The short way is to substitute the scale factors into the four general formulas of Orthogonal curvilinear coordinates: all the work was done once and for all there, and here it is arithmetic on products of $1$ , $r$ and $\sin\theta$ . The long way is to change variables in $\partial_x^2 + \partial_y^2 + \partial_z^2$ directly with the chain rule, which checks the general formula rather than assuming it. Done naively the second route is unpleasant in three dimensions, because every first derivative produces five terms and every second derivative twenty-five. Strauss's device avoids that: the spherical Laplacian is assembled out of two plane polar calculations that share the cylindrical radius $s = \sqrt{x^2+y^2}$ — one in the $(x,y)$ plane with angle $\varphi$ , one in the $(s,z)$ half-plane with angle $\theta$ — so that only the two-dimensional computation ever has to be done, and it has to be done only once.

Proof(First derivation: from the general curvilinear formulas)AdaptedAdapted from Riley, Hobson & Bence (2006), §10.10, Table 10.4, p. 369; Arfken, Weber & Harris (2013), Eqs. (3.138)–(3.143), pp. 185–187

Both systems are orthogonal, since in each case the three tangent vectors listed in the definitions are mutually perpendicular, and both are labelled in right-handed order. So the formulas of Orthogonal curvilinear coordinates apply with the scale factors recorded there.

Cylindrical. With $h_\rho = 1$ , $h_\varphi = \rho$ , $h_z = 1$ we have $H = \rho$ , and the three coefficients appearing in the Laplacian are

\frac{h_\varphi h_z}{h_\rho} = \rho, \qquad \frac{h_z h_\rho}{h_\varphi} = \frac{1}{\rho}, \qquad \frac{h_\rho h_\varphi}{h_z} = \rho .

Substituting into the general Laplacian and dividing through by $H = \rho$ ,

\nabla^2\Phi = \frac{1}{\rho}\left[ \frac{\partial}{\partial\rho}\!\left(\rho\frac{\partial\Phi}{\partial\rho}\right) + \frac{\partial}{\partial\varphi}\!\left(\frac{1}{\rho}\frac{\partial\Phi}{\partial\varphi}\right) + \frac{\partial}{\partial z}\!\left(\rho\frac{\partial\Phi}{\partial z}\right) \right] = \frac{1}{\rho}\frac{\partial}{\partial\rho}\!\left(\rho\frac{\partial\Phi}{\partial\rho}\right) + \frac{1}{\rho^2}\frac{\partial^2\Phi}{\partial\varphi^2} + \frac{\partial^2\Phi}{\partial z^2} ,

where the second and third terms simplify because $\rho$ does not depend on $\varphi$ or $z$ . The gradient, divergence and curl follow in the same way: in the divergence, $h_\varphi h_z = \rho$ , $h_z h_\rho = 1$ and $h_\rho h_\varphi = \rho$ , and the overall $1/\rho$ then cancels the $\rho$ in the last term.

Spherical. With $h_r = 1$ , $h_\theta = r$ , $h_\varphi = r\sin\theta$ we have $H = r^2\sin\theta$ , and

\frac{h_\theta h_\varphi}{h_r} = r^2\sin\theta, \qquad \frac{h_\varphi h_r}{h_\theta} = \sin\theta, \qquad \frac{h_r h_\theta}{h_\varphi} = \frac{1}{\sin\theta} .

Dividing by $H$ and using that $r$ and $\sin\theta$ are each constant with respect to the other two coordinates,

\nabla^2\Phi = \frac{1}{r^2\sin\theta}\left[ \frac{\partial}{\partial r}\!\left(r^2\sin\theta\,\frac{\partial\Phi}{\partial r}\right) + \frac{\partial}{\partial\theta}\!\left(\sin\theta\frac{\partial\Phi}{\partial\theta}\right) + \frac{\partial}{\partial\varphi}\!\left(\frac{1}{\sin\theta}\frac{\partial\Phi}{\partial\varphi}\right) \right] ,

and pulling the constant $\sin\theta$ out of the first term and $1/\sin\theta$ out of the third gives the stated formula. The divergence likewise uses $h_\theta h_\varphi = r^2\sin\theta$ , $h_\varphi h_r = r\sin\theta$ and $h_r h_\theta = r$ .

The radial identity. Expanding both sides, $\frac{1}{r^2}(r^2\Phi_r)_r = \Phi_{rr} + \frac{2}{r}\Phi_r$ , while $\frac{1}{r}(r\Phi)_{rr} = \frac{1}{r}(\Phi + r\Phi_r)_r = \frac{1}{r}(2\Phi_r + r\Phi_{rr})$ , which is the same. $\square$

Proof(Second derivation: directly, by the chain rule)AdaptedAdapted from Strauss (2008), §6.1, pp. 156–160

Let $u$ be $C^2$ . We compute $\nabla^2 u$ without using the general curvilinear formula at all.

Step 1: the plane. Let $(\rho,\varphi)$ be plane polars, $x = \rho\cos\varphi$ , $y = \rho\sin\varphi$ . Inverting the Jacobian matrix of this map gives $\partial\rho/\partial x = \cos\varphi$ , $\partial\rho/\partial y = \sin\varphi$ , $\partial\varphi/\partial x = -\sin\varphi/\rho$ and $\partial\varphi/\partial y = \cos\varphi/\rho$ — note that $\partial\rho/\partial x$ is not the reciprocal of $\partial x/\partial\rho$ . Hence the operator identities

\frac{\partial}{\partial x} = \cos\varphi\,\frac{\partial}{\partial\rho} - \frac{\sin\varphi}{\rho}\frac{\partial}{\partial\varphi}, \qquad \frac{\partial}{\partial y} = \sin\varphi\,\frac{\partial}{\partial\rho} + \frac{\cos\varphi}{\rho}\frac{\partial}{\partial\varphi} .

Apply each one twice, remembering that it also differentiates its own $\varphi$ - and $\rho$ -dependent coefficients. Sorting the outcome by which derivative of $u$ it multiplies, the two results are

	$u_{\rho\rho}$	$u_{\rho\varphi}$	$u_{\varphi\varphi}$	$u_\rho$	$u_\varphi$
$u_{xx}$	$\cos^2\varphi$	$-\dfrac{2\sin\varphi\cos\varphi}{\rho}$	$\dfrac{\sin^2\varphi}{\rho^2}$	$\dfrac{\sin^2\varphi}{\rho}$	$\dfrac{2\sin\varphi\cos\varphi}{\rho^2}$
$u_{yy}$	$\sin^2\varphi$	$+\dfrac{2\sin\varphi\cos\varphi}{\rho}$	$\dfrac{\cos^2\varphi}{\rho^2}$	$\dfrac{\cos^2\varphi}{\rho}$	$-\dfrac{2\sin\varphi\cos\varphi}{\rho^2}$

Adding the rows, the two mixed columns cancel exactly — they are equal and opposite, which is the analytic trace of the two directions being perpendicular — and the other three columns sum to $1$ , $1/\rho^2$ and $1/\rho$ by $\cos^2 + \sin^2 = 1$ . Therefore

u_{xx} + u_{yy} = u_{\rho\rho} + \frac{1}{\rho}u_\rho + \frac{1}{\rho^2}u_{\varphi\varphi} . \tag{$*$}

Step 2: cylindrical polars. The cylindrical coordinates leave $z$ untouched, so adding $u_{zz}$ to $(*)$ gives the cylindrical Laplacian immediately, and $u_{\rho\rho} + \rho^{-1}u_\rho = \rho^{-1}(\rho u_\rho)_\rho$ puts it in the tabulated form.

Step 3: spherical polars, by using $(*)$ twice. Write $s = \sqrt{x^2+y^2}$ for the distance from the $z$ -axis, so that $u$ may be regarded as a function of $(s,\varphi,z)$ or of $(r,\theta,\varphi)$ , with

s = r\sin\theta, \qquad z = r\cos\theta, \qquad r = \sqrt{s^2+z^2} .

Two plane polar systems are now visible.

In the half-plane $\varphi = \text{const}$ , the pair $(r,\theta)$ is plane polar for the Cartesian-like pair $(z,s)$ , with $\theta$ measured from the $z$ -axis. Applying $(*)$ with $z$ in the role of $x$ and $s$ in the role of $y$ : $u_{zz} + u_{ss} = u_{rr} + \frac{1}{r}u_r + \frac{1}{r^2}u_{\theta\theta} .$
In the plane $z = \text{const}$ , the pair $(s,\varphi)$ is plane polar for $(x,y)$ . Applying $(*)$ again: $u_{xx} + u_{yy} = u_{ss} + \frac{1}{s}u_s + \frac{1}{s^2}u_{\varphi\varphi} .$

In both lines $u_{ss}$ means $\partial^2u/\partial s^2$ with $\varphi$ and $z$ held fixed, so it is the same quantity and may be cancelled when the two lines are added. Adding and subtracting one copy of $u_{ss}$ ,

\nabla^2 u = u_{xx}+u_{yy}+u_{zz} = u_{rr} + \frac{1}{r}u_r + \frac{1}{r^2}u_{\theta\theta} + \frac{1}{s}u_s + \frac{1}{s^2}u_{\varphi\varphi} .

Step 4: eliminate $s$ . Only $u_s$ and the explicit $s$ 's remain. Since $r = \sqrt{s^2+z^2}$ and $\theta = \arctan(s/z)$ at fixed $z$ , we have $\partial r/\partial s = s/r = \sin\theta$ and $\partial\theta/\partial s = z/(s^2+z^2) = \cos\theta/r$ , while $\varphi$ does not change with $s$ . So the chain rule gives $u_s = \sin\theta\,u_r + (\cos\theta/r)\,u_\theta$ and

\frac{1}{s}u_s = \frac{1}{r\sin\theta}\left( \sin\theta\,u_r + \frac{\cos\theta}{r}u_\theta \right) = \frac{1}{r}u_r + \frac{\cos\theta}{r^2\sin\theta}u_\theta , \qquad \frac{1}{s^2}u_{\varphi\varphi} = \frac{1}{r^2\sin^2\theta}u_{\varphi\varphi} .

Substituting,

\nabla^2u = u_{rr} + \frac{2}{r}u_r + \frac{1}{r^2}\left( u_{\theta\theta} + \cot\theta\,u_\theta \right) + \frac{1}{r^2\sin^2\theta}u_{\varphi\varphi} .

Finally $u_{rr} + \frac{2}{r}u_r = \frac{1}{r^2}(r^2u_r)_r$ and $u_{\theta\theta} + \cot\theta\,u_\theta = \frac{1}{\sin\theta}(\sin\theta\,u_\theta)_\theta$ , which is the tabulated spherical Laplacian. It agrees with the first derivation. $\square$

RemarkOriginal

Conventions, and where the formulas stop working.

The names of the angles. We follow Riley, Hobson and Bence and Arfken: $\theta$ polar, $\varphi$ azimuthal, so that $h_\varphi = r\sin\theta$ and $dV = r^2\sin\theta\,dr\,d\theta\,d\varphi$ . Many mathematics texts exchange the two letters, and Strauss warns of this explicitly on p. 158 while himself using the same convention as ours. The formulas are unchanged in substance, but a table copied from a source with the other convention and used with this one swaps $\sin\theta$ for $\sin\varphi$ and is simply wrong. Always check which angle the factor $\sin$ attaches to.
The radial symbol. $\rho$ is the cylindrical radius here and $r$ the spherical one; they differ by $\rho = r\sin\theta$ . Sources that write $r$ for both (common in two-dimensional work) invite an error when a plane result is lifted into three dimensions: the plane polar Laplacian has $\rho^{-1}(\rho\Phi_\rho)_\rho$ , the spherical one $r^{-2}(r^2\Phi_r)_r$ , and the difference in the exponent is the whole difference between $\log\rho$ and $1/r$ as the radial harmonic function.
Degenerate sets. Each table holds only where every scale factor is positive: off the $z$ -axis for cylindrical, off the origin and off the polar axis for spherical. A field perfectly smooth at the origin can have components that are undefined there — $\hat{\mathbf{e}}_r$ itself has no limit at $r = 0$ — so a singularity in a component is not by itself a singularity in the field. Conversely a formula verified only for $r > 0$ says nothing about the origin: $\nabla^2(1/r) = 0$ holds for every $r > 0$ , yet $1/r$ is not harmonic on $\mathbb{R}^3$ , since as a distribution $\nabla^2(1/r) = -4\pi\delta$ (The Dirac delta in calculations).
Why the vector Laplacian is not componentwise. $\nabla^2\mathbf{a}$ is defined by $\nabla(\nabla\cdot\mathbf{a}) - \nabla\times(\nabla\times\mathbf{a})$ , and in Cartesian coordinates this reduces to applying the scalar operator to each component, because the basis vectors are constant. Riley, Hobson and Bence make the point (p. 357) that cylindrical and spherical polars fall in the other class: the derivatives of the rotating basis vectors contribute, and the extra terms in Arfken's Eqs. (3.151) and (3.160) are exactly those contributions. The one exception is $(\nabla^2\mathbf{a})_z$ in cylindrical coordinates, because $\hat{\mathbf{e}}_z$ is constant — and that exception is what makes the pipe-flow example below so short.

Example(Laminar flow in a circular pipe)Original

Water ( $\mu = 1.0\times10^{-3}\ \mathrm{Pa\,s}$ , $\rho_{\mathrm{m}} = 998\ \mathrm{kg\,m^{-3}}$ ) flows steadily down a straight pipe of radius $a = 5.0\ \mathrm{mm}$ along the $z$ -axis. Take the velocity field

\mathbf{v} = v(\rho)\,\hat{\mathbf{e}}_z, \qquad v(\rho) = V\left(1 - \frac{\rho^2}{a^2}\right), \qquad V = 0.20\ \mathrm{m\,s^{-1}} ,

which vanishes at the wall as a viscous fluid must. Everything below uses the cylindrical table.

Incompressibility. Only $a_z = v(\rho)$ is non-zero, and it does not depend on $z$ , so every term of the divergence vanishes: $\nabla\cdot\mathbf{v} = \frac{1}{\rho}\partial_\rho(\rho\cdot 0) + \frac{1}{\rho}\partial_\varphi(0) + \partial_z v(\rho) = 0$ . The field is admissible for an incompressible fluid.
Vorticity. In the curl determinant only the $\rho\hat{\mathbf{e}}_\varphi$ cofactor survives, giving $\nabla\times\mathbf{v} = -\dfrac{dv}{d\rho}\,\hat{\mathbf{e}}_\varphi = \dfrac{2V\rho}{a^2}\,\hat{\mathbf{e}}_\varphi$ . It vanishes on the axis and is largest at the wall, where its magnitude is $2V/a = 2(0.20)/(0.0050) = 80\ \mathrm{s^{-1}}$ .
Viscous term. Because $\mathbf{v}$ points along the constant direction $\hat{\mathbf{e}}_z$ , the remark's exception applies and $(\nabla^2\mathbf{v})_z = \nabla^2 v$ , the scalar Laplacian. With no $\varphi$ - or $z$ -dependence, $\nabla^2 v = \frac{1}{\rho}\frac{d}{d\rho}\!\left(\rho\frac{dv}{d\rho}\right) = \frac{1}{\rho}\frac{d}{d\rho}\!\left(-\frac{2V\rho^2}{a^2}\right) = -\frac{4V}{a^2} ,$ a constant.
Pressure gradient. The convective term $(\mathbf{v}\cdot\nabla)\mathbf{v} = v\,\partial_z(v\hat{\mathbf{e}}_z)$ vanishes, since neither $v$ nor $\hat{\mathbf{e}}_z$ depends on $z$ . Steady Navier–Stokes therefore reduces to $\nabla p = \mu\nabla^2\mathbf{v}$ , whose only non-trivial component is $\frac{dp}{dz} = -\frac{4\mu V}{a^2} = -\frac{4(1.0\times10^{-3})(0.20)}{(5.0\times10^{-3})^2} = -32\ \mathrm{Pa\,m^{-1}} .$ The pressure falls by $32$ Pa per metre of pipe.
Flow rate. Using the area element $\rho\,d\rho\,d\varphi$ on the face $z = \text{const}$ , $Q = \int_0^{2\pi}\!\!\int_0^a V\left(1-\frac{\rho^2}{a^2}\right)\rho\,d\rho\,d\varphi = 2\pi V\left(\frac{a^2}{2} - \frac{a^2}{4}\right) = \frac{\pi V a^2}{2} = 7.85\times10^{-6}\ \mathrm{m^3\,s^{-1}} ,$ that is $7.85\ \mathrm{mL\,s^{-1}}$ , and the mean speed is $Q/(\pi a^2) = V/2 = 0.10\ \mathrm{m\,s^{-1}}$ .
Independent check. The Hagen–Poiseuille law $Q = \pi a^4\lvert dp/dz\rvert/(8\mu)$ gives $\lvert dp/dz\rvert = 8\mu Q/(\pi a^4) = 8(1.0\times10^{-3})(7.854\times10^{-6})/(\pi(5.0\times10^{-3})^4) = 32.0\ \mathrm{Pa\,m^{-1}}$ , which matches step 4. The Reynolds number $\rho_{\mathrm{m}}(V/2)(2a)/\mu = 998$ is well below the transition near $2300$ , so the laminar profile assumed at the start is self-consistent.

Example(A conducting sphere in a uniform field)Original

An uncharged conducting sphere of radius $a = 2.0\ \mathrm{cm}$ sits in a uniform electrostatic field $\mathbf{E}_\infty = E_0\,\mathbf{k}$ with $E_0 = 1.00\times10^{5}\ \mathrm{V\,m^{-1}}$ . Take as a trial potential, for $r \ge a$ ,

\Phi(r,\theta) = -E_0\left( r - \frac{a^3}{r^2} \right)\cos\theta .

It is harmonic. Write $\Phi = f(r)\cos\theta$ with $f(r) = -E_0(r - a^3r^{-2})$ . The angular term of the spherical Laplacian needs $\frac{1}{\sin\theta}\frac{d}{d\theta}\!\left(\sin\theta\,\frac{d}{d\theta}\cos\theta\right) = \frac{1}{\sin\theta}\frac{d}{d\theta}\bigl(-\sin^2\theta\bigr) = -2\cos\theta ,$ and there is no $\varphi$ -dependence, so $\nabla^2\Phi = \bigl[\frac{1}{r^2}(r^2f')' - \frac{2}{r^2}f\bigr]\cos\theta$ . Now $f' = -E_0(1 + 2a^3r^{-3})$ , so $r^2f' = -E_0(r^2 + 2a^3r^{-1})$ and $(r^2f')' = -E_0(2r - 2a^3r^{-2})$ . Dividing by $r^2$ gives $-2E_0(r^{-1} - a^3r^{-4})$ , while $2f/r^2 = -2E_0(r^{-1} - a^3 r^{-4})$ as well. The two cancel: $\nabla^2\Phi = 0$ for $r > 0$ .
It satisfies the boundary conditions. At $r = a$ the bracket vanishes, so $\Phi = 0$ on the sphere, which is therefore an equipotential as a conductor must be. As $r\to\infty$ the second term dies and $\Phi \to -E_0 r\cos\theta = -E_0 z$ , the potential of the uniform field. Since the exterior Dirichlet problem has a unique bounded solution, this is the field.
The field. From the spherical gradient, $\mathbf{E} = -\nabla\Phi$ has $E_r = -\frac{\partial\Phi}{\partial r} = E_0\left(1 + \frac{2a^3}{r^3}\right)\cos\theta, \qquad E_\theta = -\frac{1}{r}\frac{\partial\Phi}{\partial\theta} = -E_0\left(1 - \frac{a^3}{r^3}\right)\sin\theta .$ At $r = a$ the tangential component $E_\theta$ vanishes identically, as it must at a conducting surface — a check on the algebra that the boundary condition in step 2 did not supply.
Surface charge. With $E_r(a,\theta) = 3E_0\cos\theta$ , the surface density is $\sigma = \varepsilon_0E_r = 3\varepsilon_0E_0\cos\theta$ , whose largest value is at the poles: $\sigma_{\max} = 3(8.854\times10^{-12})(1.00\times10^{5}) = 2.66\times10^{-6}\ \mathrm{C\,m^{-2}} = 2.66\ \mu\mathrm{C\,m^{-2}} .$ The field at the poles is three times the applied field, independently of the radius.
Two consistency checks, both using $dA = a^2\,d\Omega$ . The total induced charge is $\oint\sigma\,dA = 3\varepsilon_0E_0a^2\int_0^{2\pi}\!\!\int_0^\pi \cos\theta\,\sin\theta\,d\theta\,d\varphi = 0 ,$ as it must be for an uncharged sphere, because $\int_0^\pi\cos\theta\sin\theta\,d\theta = 0$ . The induced dipole moment along $\mathbf{k}$ is $p = \oint \sigma\,(a\cos\theta)\,dA = 3\varepsilon_0E_0a^3\,(2\pi)\int_0^\pi\cos^2\theta\sin\theta\,d\theta = 3\varepsilon_0E_0a^3(2\pi)\tfrac23 = 4\pi\varepsilon_0a^3E_0 ,$ numerically $8.90\times10^{-11}\ \mathrm{C\,m}$ . The potential of a point dipole of that strength is $p\cos\theta/(4\pi\varepsilon_0r^2) = E_0a^3\cos\theta/r^2$ , which is exactly the second term of $\Phi$ : outside the sphere, the conductor's whole effect is that of an induced dipole at its centre.

ApplicationOriginal

Separation of variables and the special functions. These two Laplacians are the starting point of most of the classical boundary-value problems of physics, because each separates. In cylindrical coordinates, substituting $\Phi = R(\rho)\,e^{im\varphi}e^{\pm kz}$ into $\nabla^2\Phi = 0$ leaves Bessel's equation for $R$ , which is why Bessel functions describe a vibrating drumhead, the modes of a circular waveguide and the temperature in a cooling cylinder. In spherical coordinates, $\Phi = R(r)Y(\theta,\varphi)$ leaves the spherical harmonic equation for $Y$ — the angular part of the Laplacian above is precisely the operator whose eigenfunctions the $Y_\ell^m$ are, with eigenvalue $-\ell(\ell+1)$ — and the Euler equation $r^2R'' + 2rR' - \ell(\ell+1)R = 0$ for $R$ , with solutions $r^\ell$ and $r^{-\ell-1}$ . The $\ell = 1$ pair $r$ and $r^{-2}$ is what the conducting-sphere example used.

Why the radial identity is worth remembering. In quantum mechanics the radial Schrödinger equation is tamed by exactly the alternative form $r^{-1}\partial_r^2(r\Phi)$ : putting $\Phi = u(r)/r$ turns the three-dimensional radial operator into the plain second derivative $u''(r)$ , so the hydrogen atom's radial problem becomes an ordinary one-dimensional eigenvalue problem with an effective potential. The same substitution reduces the spherically symmetric heat and wave equations to one-dimensional ones, which is how the spherical wave $f(r-ct)/r$ is found.

Everyday engineering. The pipe-flow example is the design calculation behind pressure drop in any laminar duct, from a capillary viscometer to a microfluidic channel; the conducting-sphere example is the model for field enhancement at a rounded electrode, and the factor $3$ at the poles is why discharge begins there. Both calculations take a few lines only because the operators had already been tabulated in coordinates matched to the geometry — which is the whole reason for this unit.