• 358,603

# Posts Tagged ‘jacobian’

## A final pre-exam update of my notes compilation for ‘PHY452H1S Basic Statistical Mechanics’, Taught by Prof. Arun Paramekanti

Posted by peeterjoot on April 22, 2013

Here’s my third update of my notes compilation for this course, including all of the following:

April 21, 2013 Fermi function expansion for thermodynamic quantities

April 20, 2013 Relativistic Fermi Gas

April 10, 2013 Non integral binomial coefficient

April 10, 2013 energy distribution around mean energy

April 09, 2013 Velocity volume element to momentum volume element

April 04, 2013 Phonon modes

April 03, 2013 BEC and phonons

April 03, 2013 Max entropy, fugacity, and Fermi gas

April 02, 2013 Bosons

April 02, 2013 Relativisitic density of states

March 28, 2013 Bosons

plus everything detailed in the description of my previous update and before.

## Velocity volume element to momentum volume element

Posted by peeterjoot on April 9, 2013

[Click here for a PDF of this post with nicer formatting (especially if my latex to wordpress script has left FORMULA DOES NOT PARSE errors.)]

# Motivation

One of the problems I attempted had integrals over velocity space with volume element $d^3\mathbf{u}$. Initially I thought that I’d need a change of variables to momentum space, and calculated the corresponding momentum space volume element. Here’s that calculation.

# Guts

We are working with a Hamiltonian

\begin{aligned}\epsilon = \sqrt{ (p c)^2 + \epsilon_0^2 },\end{aligned} \hspace{\stretch{1}}(1.1)

where the rest energy is

\begin{aligned}\epsilon_0 = m c^2.\end{aligned} \hspace{\stretch{1}}(1.2)

Hamilton’s equations give us

\begin{aligned}u_\alpha = \frac{ p_\alpha/c^2 }{\epsilon},\end{aligned} \hspace{\stretch{1}}(1.3)

or

\begin{aligned}p_\alpha = \frac{ m u_\alpha }{\sqrt{1 - \mathbf{u}^2/c^2}}.\end{aligned} \hspace{\stretch{1}}(1.4)

This is enough to calculate the Jacobian for our volume element change of variables

\begin{aligned}du_x \wedge du_y \wedge du_z &= \frac{\partial(u_x, u_y, u_z)}{\partial(p_x, p_y, p_z)}dp_x \wedge dp_y \wedge dp_z \\ &= \frac{1}{{c^6 \left( { m^2 + (\mathbf{p}/c)^2 } \right)^{9/2}}}\begin{vmatrix}m^2 c^2 + p_y^2 + p_z^2 & - p_y p_x & - p_z p_x \\ -p_x p_y & m^2 c^2 + p_x^2 + p_z^2 & - p_z p_y \\ -p_x p_z & -p_y p_z & m^2 c^2 + p_x^2 + p_y^2\end{vmatrix}dp_x \wedge dp_y \wedge dp_z \\ &= m^2 \left( { m^2 + \mathbf{p}^2/c^2 } \right)^{-5/2}dp_x \wedge dp_y \wedge dp_z.\end{aligned} \hspace{\stretch{1}}(1.5)

That final simplification of the determinant was a little hairy, but yielded nicely to Mathematica.

Our final result for the velocity volume element in momentum space, in terms of the particle energy is

\begin{aligned}d^3 \mathbf{u} = \frac{c^6 \epsilon_0^2 } {\epsilon^5} d^3 \mathbf{p}.\end{aligned} \hspace{\stretch{1}}(1.6)

## An updated compilation of notes, for ‘PHY452H1S Basic Statistical Mechanics’, Taught by Prof. Arun Paramekanti

Posted by peeterjoot on March 3, 2013

That compilation now all of the following too (no further updates will be made to any of these) :

February 28, 2013 Rotation of diatomic molecules

February 28, 2013 Helmholtz free energy

February 26, 2013 Statistical and thermodynamic connection

February 24, 2013 Ideal gas

February 16, 2013 One dimensional well problem from Pathria chapter II

February 15, 2013 1D pendulum problem in phase space

February 14, 2013 Continuing review of thermodynamics

February 13, 2013 Lightning review of thermodynamics

February 11, 2013 Cartesian to spherical change of variables in 3d phase space

February 10, 2013 n SHO particle phase space volume

February 10, 2013 Change of variables in 2d phase space

February 10, 2013 Some problems from Kittel chapter 3

February 07, 2013 Midterm review, thermodynamics

February 06, 2013 Limit of unfair coin distribution, the hard way

February 05, 2013 Ideal gas and SHO phase space volume calculations

February 03, 2013 One dimensional random walk

February 02, 2013 1D SHO phase space

February 02, 2013 Application of the central limit theorem to a product of random vars

January 31, 2013 Liouville’s theorem questions on density and current

January 30, 2013 State counting

## Cartesian to spherical change of variables in 3d phase space

Posted by peeterjoot on February 11, 2013

[Click here for a PDF of this post with nicer formatting]

## Question: Cartesian to spherical change of variables in 3d phase space

[1] problem 2.2 (a). Try a spherical change of vars to verify explicitly that phase space volume is preserved.

## Answer

Our kinetic Lagrangian in spherical coordinates is

\begin{aligned}\mathcal{L} &= \frac{1}{{2}} m \left( \dot{r} \hat{\mathbf{r}} + r \sin\theta \dot{\phi} \hat{\boldsymbol{\phi}} + r \dot{\theta} \hat{\boldsymbol{\theta}} \right)^2 \\ &= \frac{1}{{2}} m \left( \dot{r}^2 + r^2 \sin^2\theta \dot{\phi}^2 + r^2 \dot{\theta}^2 \right)\end{aligned} \hspace{\stretch{1}}(1.0.1)

We read off our canonical momentum

\begin{aligned}p_r &= \frac{\partial {\mathcal{L}}}{\partial {r}} \\ &= m \dot{r}\end{aligned} \hspace{\stretch{1}}(1.0.2a)

\begin{aligned}p_\theta &= \frac{\partial {\mathcal{L}}}{\partial {\theta}} \\ &= m r^2 \dot{\theta}\end{aligned} \hspace{\stretch{1}}(1.0.2b)

\begin{aligned}p_\phi &= \frac{\partial {\mathcal{L}}}{\partial {\phi}} \\ &= m r^2 \sin^2\theta \dot{\phi},\end{aligned} \hspace{\stretch{1}}(1.0.2c)

and can now express the Hamiltonian in spherical coordinates

\begin{aligned}H &= \frac{1}{{2}} m \left(\left( \frac{p_r}{m} \right)^2+ r^2 \sin^2\theta \left( \frac{p_\phi}{m r^2 \sin^2\theta} \right)+ r^2 \left( \frac{p_\theta}{m r^2} \right)\right) \\ &= \frac{p_r^2}{2m} + \frac{p_\phi^2}{2 m r^2 \sin^2\theta} + \frac{p_\theta^2}{2 m r^2}\end{aligned} \hspace{\stretch{1}}(1.0.3)

Now we want to do a change of variables. The coordinates transform as

\begin{aligned}x = r \sin\theta \cos\phi\end{aligned} \hspace{\stretch{1}}(1.0.4a)

\begin{aligned}y = r \sin\theta \sin\phi\end{aligned} \hspace{\stretch{1}}(1.0.4b)

\begin{aligned}z = r \cos\theta,\end{aligned} \hspace{\stretch{1}}(1.0.4c)

or

\begin{aligned}r = \sqrt{x^2 + y^2 + z^2}\end{aligned} \hspace{\stretch{1}}(1.0.5a)

\begin{aligned}\theta = \arccos(z/r)\end{aligned} \hspace{\stretch{1}}(1.0.5b)

\begin{aligned}\phi = \arctan(y/x).\end{aligned} \hspace{\stretch{1}}(1.0.5c)

It’s not too hard to calculate the change of variables for the momenta (verified in sphericalPhaseSpaceChangeOfVars.nb). We have

\begin{aligned}p_r = \frac{x p_x + y p_y + z p_z}{\sqrt{x^2 + y^2 + z^2}}\end{aligned} \hspace{\stretch{1}}(1.0.6a)

\begin{aligned}p_\theta = \frac{(p_x x + p_y y) z - p_z (x^2 + y^2)}{\sqrt{x^2 + y^2}}\end{aligned} \hspace{\stretch{1}}(1.0.6b)

\begin{aligned}p_\phi = x p_y - y p_x\end{aligned} \hspace{\stretch{1}}(1.0.6c)

Now let’s compute the volume element in spherical coordinates. This is

\begin{aligned}d\omega &= dr d\theta d\phi dp_r dp_\theta dp_\phi \\ &= \frac{\partial(r, \theta, \phi, p_r, p_\theta, p_\phi)}{\partial(x, y, z, p_x, p_y, p_z)}dx dy dz dp_x dp_y dp_z \\ &= \begin{vmatrix} \frac{x}{\sqrt{x^2+y^2+z^2}} & \frac{y}{\sqrt{x^2+y^2+z^2}} & \frac{z}{\sqrt{x^2+y^2+z^2}} & 0 & 0 & 0 \\ \frac{x z}{\sqrt{x^2+y^2} \left(x^2+y^2+z^2\right)} & \frac{y z}{\sqrt{x^2+y^2} \left(x^2+y^2+z^2\right)} & -\frac{\sqrt{x^2+y^2}}{x^2+y^2+z^2} & 0 & 0 & 0 \\ -\frac{y}{x^2+y^2} & \frac{x}{x^2+y^2} & 0 & 0 & 0 & 0 \\ \frac{\left(y^2+z^2\right) p_x-x y p_y-x z p_z}{\left(x^2+y^2+z^2\right)^{3/2}} & \frac{\left(x^2+z^2\right) p_y-y \left(x p_x+z p_z\right)}{\left(x^2+y^2+z^2\right)^{3/2}} & \frac{\left(x^2+y^2\right) p_z-z \left(x p_x+y p_y\right)}{\left(x^2+y^2+z^2\right)^{3/2}} & \frac{x}{\sqrt{x^2+y^2+z^2}} & \frac{y}{\sqrt{x^2+y^2+z^2}} & \frac{z}{\sqrt{x^2+y^2+z^2}} \\ \frac{y z \left(y p_x-x p_y\right)-x \left(x^2+y^2\right) p_z}{\left(x^2+y^2\right)^{3/2}} & \frac{x z \left(x p_y-y p_x\right)-y \left(x^2+y^2\right) p_z}{\left(x^2+y^2\right)^{3/2}} & \frac{x p_x+y p_y}{\sqrt{x^2+y^2}} & \frac{x z}{\sqrt{x^2+y^2}} & \frac{y z}{\sqrt{x^2+y^2}} & -\sqrt{x^2+y^2} \\ p_y & -p_x & 0 & -y & x & 0 \\ \end{vmatrix}dx dy dz dp_x dp_y dp_z \\ &= dx dy dz dp_x dp_y dp_z\end{aligned} \hspace{\stretch{1}}(1.0.7)

This also has a unit determinant, as we found in the similar cylindrical change of phase space variables.

# References

[1] RK Pathria. Statistical mechanics. Butterworth Heinemann, Oxford, UK, 1996.

## Change of variables in 2d phase space

Posted by peeterjoot on February 10, 2013

[Click here for a PDF of this post with nicer formatting]

# Motivation

In [1] problem 2.2, it’s suggested to try a spherical change of vars to verify explicitly that phase space volume is preserved, and to explore some related ideas. As a first step let’s try a similar, but presumably easier change of variables, going from Cartesian to cylindrical phase spaces.

# Canonical momenta and Hamiltonian

Our cylindrical velocity is

\begin{aligned}\mathbf{v} = \dot{r} \hat{\mathbf{r}} + r \dot{\theta} \hat{\boldsymbol{\theta}},\end{aligned} \hspace{\stretch{1}}(1.2.1)

so a purely kinetic Lagrangian would be

\begin{aligned}\mathcal{L} &= \frac{1}{{2}} m \mathbf{v}^2 \\ &= \frac{1}{{2}} m \left( \dot{r}^2 + r^2 \dot{\theta}^2 \right).\end{aligned} \hspace{\stretch{1}}(1.2.2)

Our canonical momenta are

\begin{subequations}

\begin{aligned}p_r &= \frac{\partial {\mathcal{L}}}{\partial {\dot{r}}} \\ &= m \dot{r}\end{aligned} \hspace{\stretch{1}}(1.0.3a)

\begin{aligned}p_\theta &= \frac{\partial {\mathcal{L}}}{\partial {\dot{\theta}}} \\ &= m r^2 \dot{\theta}.\end{aligned} \hspace{\stretch{1}}(1.0.3b)

\end{subequations}

and our kinetic energy is

\begin{aligned}H &= \mathcal{L} \\ &= \frac{1}{{2m}} p_r^2 + \frac{1}{{2 m r^2}} p_\theta^2.\end{aligned} \hspace{\stretch{1}}(1.0.4)

Now we need to express our momenta in terms of the Cartesian coordinates. We have for the radial momentum

\begin{aligned}p_r &= m \dot{r} \\ &= m \frac{d{{}}}{dt} \sqrt{x^2 + y^2} \\ &= \frac{1}{{2}} \frac{2 m}{r} \left( x \dot{x} + y \dot{y} \right)\end{aligned} \hspace{\stretch{1}}(1.0.5)

or

\begin{aligned}p_r = \frac{1}{{r}} \left( x p_x + y p_y \right)\end{aligned} \hspace{\stretch{1}}(1.0.6)

\begin{aligned}p_\theta &= m r^2 \frac{d{{\theta}}}{dt} \\ &= m r^2 \frac{d{{}}}{dt} \arctan \left( \frac{y}{x} \right).\end{aligned} \hspace{\stretch{1}}(1.0.7)

After some reduction (cyclindrialMomenta.nb), we find

\begin{aligned}p_\theta = p_y x - p_x y.\end{aligned} \hspace{\stretch{1}}(1.0.8)

We can assemble these into a complete set of change of variable equations

\begin{subequations}

\begin{aligned}r = \sqrt{x^2 + y^2}\end{aligned} \hspace{\stretch{1}}(1.0.9a)

\begin{aligned}\theta = \arctan\left( \frac{y}{x} \right)\end{aligned} \hspace{\stretch{1}}(1.0.9b)

\begin{aligned}p_r = \frac{1}{{\sqrt{x^2 + y^2}}} \left( x p_x + y p_y \right)\end{aligned} \hspace{\stretch{1}}(1.0.9c)

\begin{aligned}p_\theta = p_y x - p_x y.\end{aligned} \hspace{\stretch{1}}(1.0.9d)

\end{subequations}

Our phase space volume element change of variables is

\begin{aligned}dr d\theta dp_r dp_\theta &= \frac{\partial(r, \theta, p_r, p_\theta)}{\partial(x, y, p_x, p_y)}dx dy dp_x dp_y \\ &= \begin{vmatrix} \frac{x}{\sqrt{x^2+y^2}} & \frac{y}{\sqrt{x^2+y^2}} & 0 & 0 \\ -\frac{y}{x^2+y^2} & \frac{x}{x^2+y^2} & 0 & 0 \\ \frac{y \left(y p_x-x p_y\right)}{\left(x^2+y^2\right)^{3/2}} & \frac{x \left(x p_y-y p_x\right)}{\left(x^2+y^2\right)^{3/2}} & \frac{x}{\sqrt{x^2+y^2}} & \frac{y}{\sqrt{x^2+y^2}} \\ p_y & -p_x & -y & x \end{vmatrix}dx dy dp_x dp_y \\ &= \frac{x^2 + y^2}{\left(x^2 + y^2\right)^{3/2}}\frac{x^2 + y^2}{\left(x^2 + y^2\right)^{1/2}} \\ &= dx dy dp_x dp_y.\end{aligned} \hspace{\stretch{1}}(1.0.10)

We see explicitly that this point transformation has a unit Jacobian, preserving area.

# References

[1] RK Pathria. Statistical mechanics. Butterworth Heinemann, Oxford, UK, 1996.

## PHY452H1S Basic Statistical Mechanics. Problem Set 3: State counting

Posted by peeterjoot on February 10, 2013

[Click here for a PDF of this post with nicer formatting]

# Disclaimer

This is an ungraded set of answers to the problems posed.

## Question: Surface area of a sphere in d-dimensions

Consider a $d$-dimensional sphere with radius $R$. Derive the volume of the sphere $V_d(R)$ and its surface area $S_d(R)$ using $S_d(R) = dV_d(R)/dR$.

## Answer

Let’s start with some terminology and notation, borrowing from [2]

1. $1$-sphere : 2D circle, with “volume” $V_2$
2. $2$-sphere : 3D sphere, with volume $V_3$
3. $3$-sphere : 4D Euclidean hypersphere, with “volume” $V_4$
4. $4$-sphere : 5D Euclidean hypersphere, with “volume” $V_5$
5. $\cdots$

To calculate the volume, we require a parameterization allowing for expression of the volume element in an easy to integrate fashion. For the $1$-sphere, we can use the usual circular and spherical coordinate volume elements

\begin{aligned}V_2(R) = 4 \int_0^R r dr \int_0^{\pi/2} d\theta = 4 \frac{R^2}{2} \frac{\pi}{2} = \pi R^2\end{aligned} \hspace{\stretch{1}}(1.0.1a)

\begin{aligned}V_3(R) = 8 \int_0^R r^2 dr \int_0^{\pi/2} \sin\theta d\theta \int_0^{\pi/2} d\phi= 8 \frac{R^3}{3} (1) \frac{\pi}{2}= \frac{4}{3} \pi R^3\end{aligned} \hspace{\stretch{1}}(1.0.1b)

Here, to simplify the integration ranges, we’ve calculated the “volume” integral for just one quadrant or octet of the circle and sphere respectively, as in (Fig 1)

Fig1: Integrating over just one quadrant of circle or octet of sphere

How will we generalize this to higher dimensions? To calculate the volume elements in a systematic fashion, we can introduce a parameterization and use a Jacobian to change from Cartesian coordinates. For the $1$-volume and $2$-volume cases those parameterizations were the familiar

\begin{aligned}\begin{aligned}x_1 &= r \cos\theta \\ x_0 &= r \sin\theta\end{aligned}\end{aligned} \hspace{\stretch{1}}(1.0.2a)

\begin{aligned}\begin{aligned}x_2 &= r \cos\theta \\ x_1 &= r \sin\theta \cos\phi \\ x_0 &= r \sin\theta \sin\phi\end{aligned}\end{aligned} \hspace{\stretch{1}}(1.0.2b)

Some reflection shows that this generalizes nicely. Let’s use shorthand

\begin{aligned}C_k = \cos\theta_k\end{aligned} \hspace{\stretch{1}}(1.0.3a)

\begin{aligned}S_k = \sin\theta_k,\end{aligned} \hspace{\stretch{1}}(1.0.3b)

and pick $V_5$, say, a dimension bigger than the 2D or 3D cases that we can do by inspection, we can parameterize with

\begin{aligned}\begin{aligned}x_4 &= r C_4 \\ x_3 &= r S_4 C_3 \\ x_2 &= r S_4 S_3 C_2 \\ x_1 &= r S_4 S_3 S_2 C_1 \\ x_0 &= r S_4 S_3 S_2 S_1.\end{aligned}\end{aligned} \hspace{\stretch{1}}(1.0.4)

Our volume element

\begin{aligned}dV_5 = \underbrace{\frac{\partial(x_4, x_3, x_2, x_1, x_0)}{\partial(r, \theta_4, \theta_3, \theta_2, \theta_1)} }_{\equiv J_5}dr d\theta_4 d\theta_3 d\theta_2 d\theta_1\end{aligned} \hspace{\stretch{1}}(1.0.5)

That Jacobian is

\begin{aligned}J_5 = \left\lvert\begin{array}{lllll}C_4 & - r S_4 & \phantom{-} 0 & \phantom{-} 0 & \phantom{-} 0 \\ S_4 C_3 & \phantom{-} r C_4 C_3 & - r S_4 S_3 & \phantom{-} 0 & \phantom{-} 0 \\ S_4 S_3 C_2 & \phantom{-} r C_4 S_3 C_2 & \phantom{-} r S_4 C_3 C_2 & - r S_4 S_3 S_2 & \phantom{-} 0 \\ S_4 S_3 S_2 C_1 & \phantom{-} r C_4 S_3 S_2 C_1 & \phantom{-} r S_4 C_3 S_2 C_1 & \phantom{-} r S_4 S_3 C_2 C_1 & - r S_4 S_3 S_2 S_1 \\ S_4 S_3 S_2 S_1 & \phantom{-} r C_4 S_3 S_2 S_1 & \phantom{-} r S_4 C_3 S_2 S_1 & \phantom{-} r S_4 S_3 C_2 S_1 & \phantom{-} r S_4 S_3 S_2 C_1 \end{array}\right\rvert=r^4 S_4^3 S_3^2 S_2^1 S_1^0\left\lvert\begin{array}{lllll} \begin{array}{l} C_4 \end{array} \begin{array}{l} \phantom{S_3 S_2 S_1} \end{array} & \begin{array}{l} -S_4 \end{array} \begin{array}{l} \phantom{S_3 S_2 S_1} \end{array} & \begin{array}{l} \phantom{-}0 \end{array} & \begin{array}{l} \phantom{-}0 \end{array} & \begin{array}{l} \phantom{-}0 \end{array}\\ \begin{array}{l}S_4 \\ S_4 \\ S_4 \\ S_4 \end{array}\boxed{\begin{array}{l} C_3 \\ S_3 C_2 \\ S_3 S_2 C_1 \\ S_3 S_2 S_1 \end{array} }&\begin{array}{l}\phantom{-}C_4 \\ \phantom{-}C_4 \\ \phantom{-}C_4 \\ \phantom{-}C_4 \end{array}\boxed{\begin{array}{l}C_3 \\ S_3 C_2 \\ S_3 S_2 C_1 \\ S_3 S_2 S_1 \end{array}}&\begin{array}{l} - S_3 \\ \phantom{-} C_3 C_2 \\ \phantom{-} C_3 S_2 C_1 \\ \phantom{-} C_3 S_2 S_1 \end{array}&\begin{array}{l}\phantom{-} 0 \\ - S_2 \\ \phantom{-} C_2 C_1 \\ \phantom{-} C_2 S_1 \end{array}&\begin{array}{l}\phantom{-}0 \\ \phantom{-}0 \\ - S_1 \\ \phantom{-}C_1 \end{array}\end{array}\right\rvert\end{aligned} \hspace{\stretch{1}}(1.0.6)

Observe above that the cofactors of both the $1,1$ and the $1,2$ elements, when expanded along the first row, have a common factor. This allows us to work recursively

\begin{aligned}J_5 = r^4 S_4^3 S_3^2 S_2^1 S_1^0(C_4^2 - - S_4^2)\left\lvert\begin{array}{llll}C_3 & - S_3 & \phantom{-} 0 & \phantom{-} 0 \\ S_3 C_2 & \phantom{-} C_3 C_2 & - S_2 & \phantom{-} 0 \\ S_3 S_2 C_1 & \phantom{-} C_3 S_2 C_1 & \phantom{-} C_2 C_1 & - S_1 \\ S_3 S_2 S_1 & \phantom{-} C_3 S_2 S_1 & \phantom{-} C_2 S_1 & \phantom{-} C_1 \end{array}\right\rvert=r S_4^3 J_4\end{aligned} \hspace{\stretch{1}}(1.0.7)

Similarly for the 4D volume

\begin{aligned}J_4 = r^3 S_3^2 S_2^1 S_1^0\left\lvert\begin{array}{llll}C_3 & - S_3 & \phantom{-} 0 & \phantom{-} 0 \\ S_3 C_2 & \phantom{-} C_3 C_2 & - S_2 & \phantom{-} 0 \\ S_3 S_2 C_1 & \phantom{-} C_3 S_2 C_1 & \phantom{-} C_2 C_1 & - S_1 \\ S_3 S_2 S_1 & \phantom{-} C_3 S_2 S_1 & \phantom{-} C_2 S_1 & \phantom{-} C_1 \end{array}\right\rvert= r^3 S_3^2 S_2^1 S_1^0 (C_2^2 + S_2^2)\left\lvert\begin{array}{lll}\phantom{-} C_2 & - S_2 & \phantom{-} 0 \\ \phantom{-} S_2 C_1 & \phantom{-} C_2 C_1 & - S_1 \\ \phantom{-} S_2 S_1 & \phantom{-} C_2 S_1 & \phantom{-} C_1 \end{array}\right\rvert= r S_3^2 J_3\end{aligned} \hspace{\stretch{1}}(1.0.8)

and for the 3D volume

\begin{aligned}J_3 = r^2 S_2^1 S_1^0\left\lvert\begin{array}{lll}\phantom{-} C_2 & - S_2 & \phantom{-} 0 \\ \phantom{-} S_2 C_1 & \phantom{-} C_2 C_1 & - S_1 \\ \phantom{-} S_2 S_1 & \phantom{-} C_2 S_1 & \phantom{-} C_1 \end{array}\right\rvert= r^2 S_2^1 S_1^0 (C_2^2 + S_2^2)\left\lvert\begin{array}{ll}C_1 & - S_1 \\ S_1 & \phantom{-} C_1 \end{array}\right\rvert= r S_2^1 J_2,\end{aligned} \hspace{\stretch{1}}(1.0.9)

and finally for the 2D volume

\begin{aligned}J_2= r S_1^0\left\lvert\begin{array}{ll}C_1 & - S_1 \\ S_1 & \phantom{-} C_1 \end{array}\right\rvert= r S_1^0.\end{aligned} \hspace{\stretch{1}}(1.0.10)

Putting all the bits together, the “volume” element in the n-D space is

\begin{aligned}dV_n = dr d\theta_{n-1} \cdots d\theta_{1} r^{n-1} \prod_{k=0}^{n-2} (\sin\theta_{k+1})^{k}, \end{aligned} \hspace{\stretch{1}}(1.0.11)

and the total volume is

\begin{aligned}V_n(R) = 2^n \int_0^R r^{n-1} dr \prod_{k=0}^{n-2} \int_0^{\pi/2} d\theta \sin^k \theta= 2^n \frac{R^{n}}{n}\prod_{k=0}^{n-2} \int_0^{\pi/2} d\theta \sin^k \theta= \frac{4 R^2}{n} (n-2) 2^{n-2} \frac{R^{n-2}}{n-2}\prod_{k=n-3}^{n-2} \int_0^{\pi/2} d\theta \sin^k \theta\prod_{k=0}^{n-4} \int_0^{\pi/2} d\theta \sin^k \theta=4 R^2 \frac{n-2}{n} V_{n-2}(R)\int_0^{\pi/2} d\theta \sin^{n-2} \theta\int_0^{\pi/2} d\theta \sin^{n-3} \theta=4 R^2 \frac{n-2}{n} V_{n-2}(R)\frac{\pi}{2 (n-2)}.\end{aligned} \hspace{\stretch{1}}(1.0.12)

Note that a single of these sine power integrals has a messy result (see: statMechProblemSet3.nb)

\begin{aligned}\int_0^{\pi/2} d\theta \sin^{k} \theta=\frac{\sqrt{\pi } \Gamma \left(\frac{k+1}{2}\right)}{2 \Gamma \left(\frac{k}{2}+1\right)},\end{aligned} \hspace{\stretch{1}}(1.0.13)

but the product is simple, and that result, computed in statMechProblemSet3.nb, is inserted above, providing an expression for the n-D volume element as a recurrence relation

\begin{aligned}\boxed{V_{n} = \frac{2 \pi}{n} R^2 V_{n-2}(R)}\end{aligned} \hspace{\stretch{1}}(1.0.14)

With this recurrence relation we can find the volume $V_{2n}$ or $V_{2n+1}$ in terms of $V_2$ and $V_3$ respectively. For the even powers we have

\begin{aligned}\begin{array}{lll}V_{2(2)} &= \frac{2 \pi R^2}{4} V_{2} &= \frac{2^1 \pi R^2}{2(2)} V_{2} \\ V_{2(3)} &= \frac{(2 \pi R^2)^2}{(6)(4)} V_{2} &= \frac{2^2 (\pi R^2)^2}{2^2 (3)(2)} V_{2} \\ V_{2(4)} &= \frac{(2 \pi R^2)^3}{(8)(6)(4)} V_{2} &= \frac{2^3 (\pi R^2)^3}{2^3 (4)(3)(2)} V_{2} \end{array}\end{aligned} \hspace{\stretch{1}}(1.0.15)

\begin{aligned}V_{2(n)} = \frac{2^n (\pi R^2)^{n-1}}{(2n)!!} V_{2} = \frac{(\pi R^2)^{n-1}}{(n)!} V_{2} = \frac{(\pi R^2)^{n-1}}{(n)!} \pi R^2\end{aligned} \hspace{\stretch{1}}(1.0.16)

This gives us a closed form expression for the even powers for $n \ge 2$

\begin{aligned}\boxed{V_{2(n)} = \frac{(\pi R^2)^{n}}{n!}.}\end{aligned} \hspace{\stretch{1}}(1.0.17)

Observe that this also conveniently gives $V_2 = \pi R^2$, so is actually valid for $n \ge 1$. Now for the odd powers

\begin{aligned}\begin{array}{lll}V_{2(2)+1} &= \frac{2 \pi R^2}{5} V_{3} &= 3 \frac{2^1 \pi R^2}{5!!} V_{3} \\ V_{2(3)+1} &= 3 \frac{(2 \pi R^2)^2}{7 5 3} V_{3} &= 3 \frac{2^2 (\pi R^2)^2}{7!!} V_{2} \\ V_{2(4)+1} &= 3 \frac{(2 \pi R^2)^3}{(9)(7)(5)(3)} V_{2} &= 3 \frac{2^3 (\pi R^2)^3}{9!!} V_{2} \end{array}\end{aligned} \hspace{\stretch{1}}(1.0.18)

\begin{aligned}V_{2(n)+1} = 3 \frac{2^{n-1} (\pi R^2)^{n-1}}{(2n +1)!!} V_{2} = 3 \frac{2^{n-1} (\pi R^2)^{n-1}}{(2n +1)!!} \frac{4}{3} \pi R^3\end{aligned} \hspace{\stretch{1}}(1.0.19)

So, for $n \ge 2$ we have

\begin{aligned}\boxed{V_{2(n)+1} = \frac{2^{n+1} \pi^n R^{2 n + 1}}{(2n +1)!!}.}\end{aligned} \hspace{\stretch{1}}(1.0.20)

As with the even powered expression 1.0.17 we see that this is also good for $n = 1$, yielding $V_3 = 4 \pi R^2/3$ as required for the 3D spherical volume.

The even and odd power expressions don’t look quite different on the surface, but can be put into a consistent form, when expressed in terms of the gamma function.

For the even powers, using a substitution $2 n = m$ we have for even values of $m \ge 2$

\begin{aligned}V_{m} = \frac{\pi^{m/2} R^m}{\Gamma(n + 1)}= \frac{\pi^{m/2} R^m}{\Gamma(m/2 + 1)}.\end{aligned} \hspace{\stretch{1}}(1.0.21)

For the even powers, with the help of [1] we find

\begin{aligned}(2 n + 1)!! = \frac{2^{n+1} \Gamma\left(n + \frac{3}{2}\right)}{\sqrt{\pi}}.\end{aligned} \hspace{\stretch{1}}(1.0.22)

This gives us

\begin{aligned}V_{2(n)+1} = \frac{\pi^{n + 1/2} R^{2 n + 1}}{\Gamma\left(n + \frac{3}{2}\right)}= \frac{\pi^{n + 1/2} R^{2 n + 1}}{\Gamma\left(\frac{(2 n + 1) + 2}{2}\right)}.\end{aligned} \hspace{\stretch{1}}(1.0.23)

Writing $m = 2 n + 1$, or $n = (m - 1)/2$ we have for odd values of $m \ge 3$ a match with the even powered expression of 1.0.21

\begin{aligned}\boxed{V_{m} = \frac{ \pi^{m/2} R^{m} }{ \Gamma\left( m/2 + 1 \right)}.}\end{aligned} \hspace{\stretch{1}}(1.0.24)

We’ve shown that this is valid for any dimension $m \ge 2$.

Tabulating some values of these for $n \in [2, 7]$ we have respectively

\begin{aligned}\pi r^2,\frac{4 \pi r^3}{3},\frac{\pi ^2 r^4}{2},\frac{8 \pi ^2 r^5}{15},\frac{\pi ^3 r^6}{6},\frac{16 \pi ^3 r^7}{105},\end{aligned} \hspace{\stretch{1}}(1.0.25)

The only task left is computation of the surface area. That comes by inspection and is

\begin{aligned}S_{m}(R) =\frac{d V_m(R)}{d R}= \frac{\pi^{m/2} m R^{m-1}}{\Gamma\left( m/2 + 1 \right)}.\end{aligned} \hspace{\stretch{1}}(1.0.26)

Again for $m \in [2, 7]$ we have

\begin{aligned}2 \pi r,4 \pi r^2,2 \pi ^2 r^3,\frac{8 \pi ^2 r^4}{3},\pi ^3 r^5,\frac{16 \pi ^3 r^6}{15}.\end{aligned} \hspace{\stretch{1}}(1.0.27)

## Question: State counting – polymer

A typical protein is a long chain molecule made of very many elementary units called amino acids – it is an example of a class of such macromolecules called polymers. Consider a protein made $N$ amino acids, and assume each amino acid is like a sphere of radius a. As a toy model assume that the protein configuration is like a random walk with each amino acid being one “step”, i.e., the center-to-center vector from one amino acid to the next is a random vector of length $2a$ and ignore any issues with overlapping spheres (so-called “excluded volume” constraints). Estimate the spatial extent of this protein in space. Typical proteins assume a compact form in order to be functional. In this case, taking the constraint of nonoverlapping spheres, estimate the radius of such a compact protein assuming it has an overall spherical shape with fully packed amino acids (ignore holes in the packing, and use only volume ratios to make this estimate). With $N = 300$ and $a \approx 5 \text{\AA}$, estimate these sizes for the random walk case as well as the compact globular case.

## Answer

We are considering a geometry like that of (Fig 2), depicted in two dimensions for ease of illustration.

Fig2: Touching “spheres”

From the geometry, if $\mathbf{c}_k$ is the vector to the center of the $k$th sphere, we have for some random unit vector $\hat{\mathbf{r}}_k$

\begin{aligned}\mathbf{c}_{k+1} = \mathbf{c}_k + 2 a \hat{\mathbf{r}}_k\end{aligned} \hspace{\stretch{1}}(1.0.28)

Proceeding recursively, writing $d_N$, and $n = N -1 > 0$, we have for the difference of the positions of the first and last centers of the chain

\begin{aligned}d_N = \left\lvert {\mathbf{c}_N - \mathbf{c}_1} \right\rvert = 2 a \left\lvert { \sum_{k = 1}^{n} \hat{\mathbf{r}}_k} \right\rvert= 2 a \left( \sum_{k,m = 1}^{n} \hat{\mathbf{r}}_k \cdot \hat{\mathbf{r}}_m \right)^{1/2}= 2 a \left( n + 2 \sum_{1 \le k < m \le n} \hat{\mathbf{r}}_k \cdot \hat{\mathbf{r}}_m \right)^{1/2}= 2 a \sqrt{n} \left( 1 + \frac{2}{n} \sum_{1 \le k < m \le n} \hat{\mathbf{r}}_k \cdot \hat{\mathbf{r}}_m \right)^{1/2}\end{aligned} \hspace{\stretch{1}}(1.0.29)

The $\hat{\mathbf{r}}_k$‘s clearly cannot be completely random since we have a constraint that $\hat{\mathbf{r}}_k \cdot \hat{\mathbf{r}}_{k+1} > -\cos\pi/3$, or else two adjacent spheres will overlap. There will also be overlapping constraints for longer portions of the chain that are harder to express. We are ignoring both such constraints, and seek the ensemble average of all systems of the form 1.0.29.

Employing random azimuthal and polar angular variables $\theta_k$, and $\phi_k$, we have

\begin{aligned}\hat{\mathbf{r}}_k = \begin{bmatrix}\sin\theta_k \cos\phi_k \\ \sin\theta_k \sin\phi_k \\ \cos\theta_k,\end{bmatrix}\end{aligned} \hspace{\stretch{1}}(1.0.30)

so that the average polymer length is

\begin{aligned}\begin{aligned}\left\langle{{d_N }}\right\rangle &=2 a \sqrt{n} \left( \frac{1}{{(2 \pi)(\pi)}} \right)^{n} \int_{\theta_j \in [0, \pi]}d\theta_1 d\theta_2 \cdots d\theta_n \int_{\phi_j \in [0, 2 \pi]}d\phi_1 d\phi_2 \cdots d\phi_n \times \\ &\quad \left( 1 + \frac{2}{n} \sum_{1 \le k < m \le n} \sin\theta_k \cos\phi_k \sin\theta_m \cos\phi_m + \sin\theta_k \sin\phi_k \sin\theta_m \sin\phi_m + \cos\theta_k \cos\theta_m \right)^{1/2}\end{aligned}\end{aligned} \hspace{\stretch{1}}(1.0.31)

Observing that even $\int \sqrt{1 + a \cos\theta} d\theta$ is an elliptic integral, we don’t have any hope of evaluating this in closed form. However, to first order, we have

\begin{aligned}\begin{aligned}d_N&\approx2 a \sqrt{n} \left( \frac{1}{{(2 \pi)(\pi)}} \right)^{n} \int_{\theta_j \in [0, \pi]}d\theta_1 d\theta_2 \cdots d\theta_n \int_{\phi_j \in [0, 2 \pi]}d\phi_1 d\phi_2 \cdots d\phi_n \times \\ &\quad \left( 1 + \frac{1}{n} \sum_{1 \le k < m \le n} \sin\theta_k \cos\phi_k \sin\theta_m \cos\phi_m + \sin\theta_k \sin\phi_k \sin\theta_m \sin\phi_m + \cos\theta_k \cos\theta_m \right) \\ &=2 a \sqrt{n} \left( \frac{1}{{(2 \pi)(\pi)}} \right)^{n} \left( 2 \pi^2 \right)^n \\ &\quad +2 a \sqrt{n} \left( \frac{1}{{(2 \pi)(\pi)}} \right)^{n} \left( 2 \pi^2 \right)^{n - 2} \left( \frac{1}{{2}} n(n+1) \right) \frac{1}{{n}} \times \\ &\quad \int_0^{\pi}d\theta \int_0^{\pi}d\theta'\int_0^{2 \pi}d\phi\int_0^{2 \pi}d\phi'\left( \sin\theta \cos\phi \sin\theta' \cos\phi' + \sin\theta \sin\phi \sin\theta' \sin\phi' + \cos\theta \cos\theta' \right)\end{aligned}\end{aligned} \hspace{\stretch{1}}(1.0.32)

The $\int_0^{2 \pi} \cos\phi$ integrals kill off the first term, the $\int_0^{2 \pi} \sin\phi$ integrals kill of the second, and the $\int_0^\pi \cos\theta$ integral kill of the last term, and we are left with just

\begin{aligned}d_N \approx 2 a \sqrt{N-1}.\end{aligned} \hspace{\stretch{1}}(1.0.33)

Ignoring the extra $2 \times a$ of the end points, and assuming that for large $N$ we have $\sqrt{N-1} \approx \sqrt{N}$, the spatial extent of the polymer chain is

\begin{aligned}\boxed{d_N \approx 2 a \sqrt{N}.}\end{aligned} \hspace{\stretch{1}}(1.0.34)

Spherical packing

Assuming the densest possible spherical packing, a face centered cubic [3], as in (Fig 3), we see that the density of such a spherical packing is

Fig3: Element of a face centered cubic

\begin{aligned}\eta_{\mathrm{FCC}} =\frac{\left( 8 \times \frac{1}{{8}} + 6 \times \frac{1}{{2}} \right) \frac{4}{3} \pi r^3}{ \left( \sqrt{8 r^2} \right)^3 }= \frac{\pi}{\sqrt{18}}.\end{aligned} \hspace{\stretch{1}}(1.0.35)

With a globular radius of $R$ and an atomic radius of $a$, and density $\eta$ we have

\begin{aligned}\eta \frac{4}{3} \pi R^3 = N \frac{4}{3} \pi a^3,\end{aligned} \hspace{\stretch{1}}(1.0.36)

so that the globular radius $R$ is

\begin{aligned}\boxed{R_N(\eta) = a \sqrt[3]{\frac{N}{\eta}}.}\end{aligned} \hspace{\stretch{1}}(1.0.37)

Some numbers

With $N = 300$ and $a \approx 5 \text{\AA}$, and ignoring spaces (i.e. $\eta = 1$, for a non-physical infinite packing), our globular diameter is approximately

\begin{aligned}2 \times 5 \text{\AA} \sqrt[3]{300} \approx 67 \text{\AA}.\end{aligned} \hspace{\stretch{1}}(1.0.38)

This is actually not much different than the maximum spherical packing of an FCC lattice, which results a slightly larger globular cluster diameter

\begin{aligned}2 \times 5 \text{\AA} \sqrt[3]{300 \sqrt{18}/\pi} \approx 74 \text{\AA}.\end{aligned} \hspace{\stretch{1}}(1.0.39)

Both however, are much less than the end to end length of the random walk polymer chain

\begin{aligned}2 (5 \text{\AA}) \sqrt{300} \approx 173 \text{\AA}.\end{aligned} \hspace{\stretch{1}}(1.0.40)

## Question: State counting – spins

Consider a toy model of a magnet where the net magnetization arises from electronic spins on each atom which can point in one of only two possible directions – Up/North or Down/South. If we have a system with $N$ spins, and if the magnetization can only take on values $\pm 1$ (Up = $+1$, Down = $-1$), how many configurations are there which have a total magnetization $m$, where $m = (N_\uparrow - N_\downarrow)/N$ (note that $N_\uparrow + N_\downarrow = N$)? Simplify this result assuming $N \gg 1$ and a generic $m$ (assume we are not interested in the extreme case of a fully magnetized system where $m = \pm 1$). Find the value of the magnetization $m$ for which the number of such microscopic states is a maximum. For $N = 20$, make a numerical plot of the number of states as a function of the magnetization (note: $-1 \le m \le 1$) without making the $N \gg 1$ assumption.

## Answer

For the first couple values of $N$, lets enumerate the spin sample spaces, their magnetization.

$N = 1$

1. $\uparrow$ : $m = 1$
2. $\downarrow$, $m = -1$

$N = 2$

1. $\uparrow \uparrow$ : $m = 1$
2. $\uparrow \downarrow$ : $m = 0$
3. $\downarrow \uparrow$, $m = 0$
4. $\downarrow \downarrow$, $m = -1$

$N = 3$

1. $\uparrow \uparrow \uparrow$ : $m = 1$
2. $\uparrow \uparrow \downarrow$ : $m = 1/3$
3. $\uparrow \downarrow \uparrow$ : $m = 1/3$
4. $\uparrow \downarrow \downarrow$ : $m = -1/3$
5. $\downarrow \uparrow \uparrow$ : $m = 1/3$
6. $\downarrow \uparrow \downarrow$ : $m = -1/3$
7. $\downarrow \downarrow \uparrow$ : $m = -1/3$
8. $\downarrow \downarrow \downarrow$ : $m = 1$

The respective multiplicities for $N = 1,2,3$ are $\{1\}$, $\{1, 2, 1\}$, $\{1, 3, 3, 1\}$. It’s clear that these are just the binomial coefficients. Let’s write for the multiplicities

\begin{aligned}g(N, m) = \binom{N}{i(m)}\end{aligned} \hspace{\stretch{1}}(1.0.41)

where $i(m)$ is a function that maps from the magnetization values $m$ to the integers $[0, N]$. Assuming

\begin{aligned}i(m) = a m + b,\end{aligned} \hspace{\stretch{1}}(1.0.42)

where $i(-1) = 0$ and $i(1) = N$, we solve

\begin{aligned}\begin{aligned}a (-1) + b &= 0 \\ a (1) + b &= N,\end{aligned}\end{aligned} \hspace{\stretch{1}}(1.0.43)

so

\begin{aligned}i(m) = \frac{N}{2}(m + 1)\end{aligned} \hspace{\stretch{1}}(1.0.44)

and

\begin{aligned}g(N, m) = \binom{N}{\frac{N}{2}(1 + m)} = \frac{N!}{\left(\frac{N}{2}(1 + m)\right)!\left(\frac{N}{2}(1 - m)\right)!}\end{aligned} \hspace{\stretch{1}}(1.0.45)

From

\begin{aligned}\begin{aligned}2 m &= N_{\uparrow} - N_{\downarrow} \\ N &= N_{\uparrow} + N_{\downarrow},\end{aligned}\end{aligned} \hspace{\stretch{1}}(1.0.46)

we see that this can also be written

\begin{aligned}\boxed{g(N, m) = \frac{N!}{(N_{\uparrow})!(N_{\downarrow})!}}\end{aligned} \hspace{\stretch{1}}(1.0.47)

Simplification for large $N$

Using Stirlings approximation

\begin{aligned}\begin{aligned}\ln g(N, m) &= \ln N! - \ln N_{\downarrow}! - \ln N_{\uparrow}! \\ &\approx \left( N + \frac{1}{{2}} \right) \ln N - \not{{N}} + \not{{ \frac{1}{{2}} \ln 2 \pi }} \\ &\quad-\left( N_{\uparrow} + \frac{1}{{2}} \right) \ln N_{\uparrow} + \not{{N_{\uparrow}}} - \not{{\frac{1}{{2}} \ln 2 \pi}} \\ &\quad-\left( N_{\downarrow} + \frac{1}{{2}} \right) \ln N_{\downarrow} + \not{{N_{\downarrow}}} - \frac{1}{{2}} \ln 2 \pi \\ &=\left( N_{\uparrow} + N_{\downarrow} + \frac{1}{{2}} + \underbrace{\frac{1}{{2}}}_{\text{add}} \right) \ln N -\frac{1}{{2}} \ln 2 \pi - \underbrace{\frac{1}{{2}} \ln N}_{\text{then subtract}} \\ &\quad -\left( N_{\uparrow} + \frac{1}{{2}} \right) \ln N_{\uparrow} \quad -\left( N_{\downarrow} + \frac{1}{{2}} \right) \ln N_{\downarrow} \\ &=-\left( N_{\uparrow} + \frac{1}{{2}} \right) \ln N_{\uparrow}/N-\left( N_{\downarrow} + \frac{1}{{2}} \right) \ln N_{\downarrow}/N- \frac{1}{{2}} \ln 2 \pi N \\ &=-\left( N_{\uparrow} + \frac{1}{{2}} \right) \ln \frac{1}{{2}}( 1 + m )-\left( N_{\downarrow} + \frac{1}{{2}} \right) \ln \frac{1}{{2}} ( 1 - m )- \frac{1}{{2}} \ln 2 \pi N \\ &\approx\left( N_{\uparrow} + \frac{1}{{2}} \right) \ln 2+\left( N_{\downarrow} + \frac{1}{{2}} \right) \ln 2 - \frac{1}{{2}} \ln 2 \pi N \\ &\quad -\left( N_{\uparrow} + \frac{1}{{2}} \right) \left(m - \frac{1}{{2}} m^2 \right) \\ &\quad -\left( N_{\downarrow} + \frac{1}{{2}} \right) \left( -m - \frac{1}{{2}} m^2 \right) \\ &= (N + 1) \ln 2 - \frac{1}{{2}} \ln 2 \pi N+ m ( N_{\downarrow} - N_{\uparrow} )+ \frac{1}{{2}} m^2 ( N + 1 ) \\ &= (N + 1) \ln 2 - \frac{1}{{2}} \ln 2 \pi N- N m^2+ \frac{1}{{2}} m^2 ( N + 1 ) \\ &\approx(N + 1) \ln 2 - \frac{1}{{2}} \ln 2 \pi N- \frac{1}{{2}} m^2 N\end{aligned}\end{aligned} \hspace{\stretch{1}}(1.0.48)

This gives us for large $N$

\begin{aligned}g(N, m) \approx 2^N \sqrt{\frac{2}{ \pi N}} e^{ -m^2 N/2 }\end{aligned} \hspace{\stretch{1}}(1.0.49)

For $N = 20$ this approximation and the exact expression are plotted in (Fig 4).

Fig4: Distribution of number of configurations for N = 20 magnets as a function of magnetization

With the large scales of this extremely peaked function, no visible difference can be seen. That difference does exist, and is plotted in (Fig 5)

Fig5: N = 20 differences between the exact binomial expression and the Gaussian approximation

# References

[1] M. Abramowitz and I.A. Stegun. \emph{Handbook of mathematical functions with formulas, graphs, and mathematical tables}, volume 55. Dover publications, 1964.

[2] Wikipedia. N-sphere — Wikipedia, The Free Encyclopedia, 2013\natexlab{a}. URL http://en.wikipedia.org/w/index.php?title=N-sphere&oldid=534164100. [Online; accessed 26-January-2013].

[3] Wikipedia. Sphere packing — wikipedia, the free encyclopedia, 2013\natexlab{b}. URL http://en.wikipedia.org/w/index.php?title=Sphere_packing&oldid=535578971. [Online; accessed 31-January-2013].

## Hypersphere volume calculation the easy way

Posted by peeterjoot on January 27, 2013

[Click here for a PDF of this post with nicer formatting (especially if my latex to wordpress script has left FORMULA DOES NOT PARSE errors.)]

# Motivation

In problem set three I solved the hypersphere volume problem the hard way. I did it the way that I thought was obvious, starting with the spherical coordinate hypervolume element for an Euclidean space. That build on the volume element calculation I’d previously done for four dimensional Euclidean and Hyperbolic spaces in [1]. This time I avoided any use of Geometric Algebra and wrote out the volume element directly using a Jacobian transformation instead of a wedge product (that leads to the same Jacobian). I then proceeded to integrate the volume element, come up with a recurrence relation for the volume, solve that for even and odd powers, and find a common expression using Gamma functions that was correct for even and odd powers. It was a labourious process, but satisfying since I’d previously tried this calculation and not succeeded.

As and after I did this calculation, which I’ll post later, I figured there had to have been an easier way. In the midterm prep reading of section 5.5 of [2] I found just that method. It’s done there in about six lines, using a trick that seemed really sneaky! I was left with a frustrated feeling, wondering how on earth somebody figured out to do it that way.

After a bit of reflection, I see that the trick is a very sensible approach. I’ll outline that here for 2 and 3 dimensions to show the technique, and the reader can generalize if desired.

# The trick

Switching to spherical or circular coordinates when there is radial symmetry isn’t anything that we could describe as trickery. For example, with $r^2 = x^2 + y^2$ in a two dimensional problem or $r^2 = x^2 + y^2 + z^2$ in a three dimensional problem, it’s clear that we’d do the following respectively if we were evaluating an integral

\begin{aligned}\iint f(r) dx dy = \iint f(r) r dr d\theta = 2 \pi \int f(r) r dr\end{aligned} \hspace{\stretch{1}}(1.0.5a)

\begin{aligned}\iiint f(r) dx dy = \iiint f(r) r^2 dr \sin d\theta d\theta d\phi = 4 \pi \int f(r) r^2 dr\end{aligned} \hspace{\stretch{1}}(1.0.5b)

In fact, for $f(r) = 1$ these give us the area and volume of the circle and sphere respectively.

So, what’s the trick? The first part is the observation that the “area” of a “volume” for the circle and the sphere are found by the derivatives of the “volumes”

\begin{aligned}\frac{d}{dr} \pi r^2 = 2 \pi r\end{aligned} \hspace{\stretch{1}}(1.0.2a)

\begin{aligned}\frac{d}{dr} \frac{4 \pi}{3} r^3 = 4 \pi r^2.\end{aligned} \hspace{\stretch{1}}(1.0.2b)

I recall being suprised that this is the case back in high school calculus. When I lost marks for not having the formula for the surface area of a sphere memorized on some test, my calculus teacher told me this title tidbit.

Back then this wasn’t obvious to me, and I complained that it was true only for the perfect symmetry of a sphere. For example, the derivative of an volume of a cube doens’t give the surface area of a cube ($d x^3/dx = 3 x^2 \ne 6 x^2$).

Once we believe or assume that the surface “area” of a hypervolume is the derivative of the volume, we can proceed with the trick. That trick is to express the volume in terms of an unknown constant. For example, for the circle and sphere the generalized “volume”s are respectively

\begin{aligned}V_2 = B r^2\end{aligned} \hspace{\stretch{1}}(1.0.3a)

\begin{aligned}V_3 = C r^3\end{aligned} \hspace{\stretch{1}}(1.0.3b)

The perimeter of the circle and surface area of the sphere are then

\begin{aligned}A_2 = 2 B r\end{aligned} \hspace{\stretch{1}}(1.0.4a)

\begin{aligned}A_3 = 3 C r^2\end{aligned} \hspace{\stretch{1}}(1.0.4b)

So, if we want to calculate integrals of the form 1.0.1.1 we can write

\begin{aligned}\iint f(r) dx dy = 2 B \int f(r) r dr\end{aligned} \hspace{\stretch{1}}(1.0.5a)

\begin{aligned}\iiint f(r) dx dy = 3 C \int f(r) r^2 dr.\end{aligned} \hspace{\stretch{1}}(1.0.5b)

The essence of the trick is to do such an integral in both Cartesian coordinates to get the left hand sides, and then do the radial right hand side integrals. Comparing these provides the constants $B$ and $C$ and thus completes the “volume” formulas for the circle and sphere.

The function chosen for this integral in the text was a Gaussian exponential $f(r) = e^{-r^2/2}$, something that is strictly radial, and can be integrated over all space. For the 2D case, we’ll integrate

\begin{aligned}\iint e^{-(x^2 +y^2)/2} dx dy = (\sqrt{2 \pi})^2= 2 B \int_0^\infty dr r e^{-r^2/2}= -2 B {\left.{e^{-r^2/2}}\right\vert}_{0}^{\infty} = 2 B.\end{aligned} \hspace{\stretch{1}}(1.0.5b)

We find that $B = \pi$, so the 2D spherical “volume” (the area of the circle) is $V = \pi r^2$.

For the 3D sphere, we have

\begin{aligned}\iiint e^{-(x^2 +y^2 + z^2)/2} dx dy = (\sqrt{2 \pi})^3= 3 C \int_0^\infty dr r^2 e^{-r^2/2}= 3 C \int_0^\infty dr e^{-r^2/2}= 3 C \sqrt{2 \pi}.\end{aligned} \hspace{\stretch{1}}(1.0.5b)

So we have for the volume of the 3D sphere, $V = 4 \pi r^3/3$ as expected. The same idea can be extended to higher dimensional spheres. The text does the even values of $N$. Treating both even and odd values, I’d expect this to yield the result I calculated with the Jacobian volume element method

\begin{aligned}V_{m} = \frac{ \pi^{m/2} R^{m} }{ \Gamma\left( m/2 + 1 \right)}.\end{aligned} \hspace{\stretch{1}}(1.0.5b)

# References

[1] Peeter Joot. Exploring physics with Geometric Algebra, chapter {Spherical and hyperspherical parametrization}. URL http://sites.google.com/site/peeterjoot/math2009/gabook.pdf.

[2] S.K. Ma. Statistical Mechanics. World Scientific, 1985. ISBN 9789971966072. URL http://books.google.ca/books?id=3XyA0SYK4XIC.

## Exploring Stokes Theorem in tensor form.

Posted by peeterjoot on February 22, 2011

# Obsolete with potential errors.

This post may be in error.  I wrote this before understanding that the gradient used in Stokes Theorem must be projected onto the tangent space of the parameterized surface, as detailed in Alan MacDonald’s Vector and Geometric Calculus.

See the post ‘stokes theorem in geometric algebra‘ [PDF], where this topic has been revisited with this in mind.

# Original Post:

[Click here for a PDF of this post with nicer formatting]

## Motivation.

I’ve worked through Stokes theorem concepts a couple times on my own now. One of the first times, I was trying to formulate this in a Geometric Algebra context. I had to resort to a tensor decomposition, and pictures, before ending back in the Geometric Algebra description. Later I figured out how to do it entirely with a Geometric Algebra description, and was able to eliminate reliance on the pictures that made the path to generalization to higher dimensional spaces unclear.

It’s my expectation that if one started with a tensor description, the proof entirely in tensor form would not be difficult. This is what I’d like to try this time. To start off, I’ll temporarily use the Geometric Algebra curl expression so I know what my tensor equation starting point will be, but once that starting point is found, we can work entirely in coordinate representation. For somebody who already knows that this is the starting point, all of this initial motivation can be skipped.

# Translating the exterior derivative to a coordinate representation.

Our starting point is a curl, dotted with a volume element of the same grade, so that the result is a scalar

\begin{aligned}\int d^n x \cdot (\nabla \wedge A).\end{aligned} \hspace{\stretch{1}}(2.1)

Here $A$ is a blade of grade $n-1$, and we wedge this with the gradient for the space

\begin{aligned}\nabla \equiv e^i \partial_i = e_i \partial^i,\end{aligned} \hspace{\stretch{1}}(2.2)

where we with with a basis (not necessarily orthonormal) $\{e_i\}$, and the reciprocal frame for that basis $\{e^i\}$ defined by the relation

\begin{aligned}e^i \cdot e_j = {\delta^i}_j.\end{aligned} \hspace{\stretch{1}}(2.3)

Our coordinates in these basis sets are

\begin{aligned}x \cdot e^i & \equiv x^i \\ x \cdot e_i & \equiv x_i\end{aligned} \hspace{\stretch{1}}(2.4)

so that

\begin{aligned}x = x^i e_i = x_i e^i.\end{aligned} \hspace{\stretch{1}}(2.6)

The operator coordinates of the gradient are defined in the usual fashion

\begin{aligned}\partial_i & \equiv \frac{\partial }{\partial {x^i}} \\ \partial^i & \equiv \frac{\partial}{\partial {x_i}}\end{aligned} \hspace{\stretch{1}}(2.7)

The volume element for the subspace that we are integrating over we will define in terms of an arbitrary parametrization

\begin{aligned}x = x(\alpha_1, \alpha_2, \cdots, \alpha_n)\end{aligned} \hspace{\stretch{1}}(2.9)

The subspace can be considered spanned by the differential elements in each of the respective curves where all but the $i$th parameter are held constant.

\begin{aligned}dx_{\alpha_i}= d\alpha_i \frac{\partial x}{\partial {\alpha_i}}= d\alpha_i \frac{\partial {x^j}}{\partial {\alpha_i}} e_j.\end{aligned} \hspace{\stretch{1}}(2.10)

We assume that the integral is being performed in a subspace for which none of these differential elements in that region are linearly dependent (i.e. our Jacobean determinant must be non-zero).

The magnitude of the wedge product of all such differential elements provides the volume of the parallelogram, or parallelepiped (or higher dimensional analogue), and is

\begin{aligned}d^n x=d\alpha_1 d\alpha_2\cdots d\alpha_n\frac{\partial x}{\partial {\alpha_n}} \wedge\cdots \wedge\frac{\partial x}{\partial {\alpha_2}}\wedge\frac{\partial x}{\partial {\alpha_1}}.\end{aligned} \hspace{\stretch{1}}(2.11)

The volume element is a oriented quantity, and may be adjusted with an arbitrary sign (or equivalently an arbitrary permutation of the differential elements in the wedge product), and we’ll see that it is convenient for the translation to tensor form, to express these in reversed order.

Let’s write

\begin{aligned}d^n \alpha = d\alpha_1 d\alpha_2 \cdots d\alpha_n,\end{aligned} \hspace{\stretch{1}}(2.12)

so that our volume element in coordinate form is

\begin{aligned}d^n x = d^n \alpha\frac{\partial {x^i}}{\partial {\alpha_1}}\frac{\partial {x^j}}{\partial {\alpha_2}}\cdots \frac{\partial {x^k}}{\partial {\alpha_{n-1}}}\frac{\partial {x^l}}{\partial {\alpha_n}}( e_l \wedge e_k \wedge \cdots \wedge e_j \wedge e_i ).\end{aligned} \hspace{\stretch{1}}(2.13)

Our curl will also also be a grade $n$ blade. We write for the $n-1$ grade blade

\begin{aligned}A = A_{b c \cdots d} (e^b \wedge e^c \wedge \cdots e^d),\end{aligned} \hspace{\stretch{1}}(2.14)

where $A_{b c \cdots d}$ is antisymmetric (i.e. $A = a_1 \wedge a_2 \wedge \cdots a_{n-1}$ for a some set of vectors $a_i, i \in 1 .. n-1$).

With our gradient in coordinate form

\begin{aligned}\nabla = e^a \partial_a,\end{aligned} \hspace{\stretch{1}}(2.15)

the curl is then

\begin{aligned}\nabla \wedge A = \partial_a A_{b c \cdots d} (e^a \wedge e^b \wedge e^c \wedge \cdots e^d).\end{aligned} \hspace{\stretch{1}}(2.16)

The differential form for our integral can now be computed by expanding out the dot product. We want

\begin{aligned}( e_l \wedge e_k \wedge \cdots \wedge e_j \wedge e_i )\cdot(e^a \wedge e^b \wedge e^c \wedge \cdots e^d)=((((( e_l \wedge e_k \wedge \cdots \wedge e_j \wedge e_i ) \cdot e^a ) \cdot e^b ) \cdot e^c ) \cdot \cdots ) \cdot e^d.\end{aligned} \hspace{\stretch{1}}(2.17)

Evaluation of the interior dot products introduces the intrinsic antisymmetry required for Stokes theorem. For example, with

\begin{aligned}( e_n \wedge e_{n-1} \wedge \cdots \wedge e_2 \wedge e_1 ) \cdot e^a a & =( e_n \wedge e_{n-1} \wedge \cdots \wedge e_3 \wedge e_2 ) (e_1 \cdot e^a) \\ & -( e_n \wedge e_{n-1} \wedge \cdots \wedge e_3 \wedge e_1 ) (e_2 \cdot e^a) \\ & +( e_n \wedge e_{n-1} \wedge \cdots \wedge e_2 \wedge e_1 ) (e_3 \cdot e^a) \\ & \cdots \\ & (-1)^{n-1}( e_{n-1} \wedge e_{n-2} \wedge \cdots \wedge e_2 \wedge e_1 ) (e_n \cdot e^a)\end{aligned}

Since $e_i \cdot e^a = {\delta_i}^a$ our end result is a completely antisymmetric set of permutations of all the deltas

\begin{aligned}( e_l \wedge e_k \wedge \cdots \wedge e_j \wedge e_i )\cdot(e^a \wedge e^b \wedge e^c \wedge \cdots e^d)={\delta^{[a}}_i{\delta^b}_j\cdots {\delta^{d]}}_l,\end{aligned} \hspace{\stretch{1}}(2.18)

and the curl integral takes it’s coordinate form

\begin{aligned}\int d^n x \cdot ( \nabla \wedge A ) =\int d^n \alpha\frac{\partial {x^i}}{\partial {\alpha_1}}\frac{\partial {x^j}}{\partial {\alpha_2}}\cdots \frac{\partial {x^k}}{\partial {\alpha_{n-1}}}\frac{\partial {x^l}}{\partial {\alpha_n}}\partial_a A_{b c \cdots d}{\delta^{[a}}_i{\delta^b}_j\cdots {\delta^{d]}}_l.\end{aligned} \hspace{\stretch{1}}(2.19)

One final contraction of the paired indexes gives us our Stokes integral in its coordinate representation

\begin{aligned}\boxed{\int d^n x \cdot ( \nabla \wedge A ) =\int d^n \alpha\frac{\partial {x^{[a}}}{\partial {\alpha_1}}\frac{\partial {x^b}}{\partial {\alpha_2}}\cdots \frac{\partial {x^c}}{\partial {\alpha_{n-1}}}\frac{\partial {x^{d]}}}{\partial {\alpha_n}}\partial_a A_{b c \cdots d}}\end{aligned} \hspace{\stretch{1}}(2.20)

We now have a starting point that is free of any of the abstraction of Geometric Algebra or differential forms. We can identify the products of partials here as components of a scalar hypervolume element (possibly signed depending on the orientation of the parametrization)

\begin{aligned}d\alpha_1 d\alpha_2\cdots d\alpha_n\frac{\partial {x^{[a}}}{\partial {\alpha_1}}\frac{\partial {x^b}}{\partial {\alpha_2}}\cdots \frac{\partial {x^c}}{\partial {\alpha_{n-1}}}\frac{\partial {x^{d]}}}{\partial {\alpha_n}}\end{aligned} \hspace{\stretch{1}}(2.21)

This is also a specific computation recipe for these hypervolume components, something that may not be obvious when we allow for general metrics for the space. We are also allowing for non-orthonormal coordinate representations, and arbitrary parametrization of the subspace that we are integrating over (our integral need not have the same dimension as the underlying vector space).

Observe that when the number of parameters equals the dimension of the space, we can write out the antisymmetric term utilizing the determinant of the Jacobian matrix

\begin{aligned}\frac{\partial {x^{[a}}}{\partial {\alpha_1}}\frac{\partial {x^b}}{\partial {\alpha_2}}\cdots \frac{\partial {x^c}}{\partial {\alpha_{n-1}}}\frac{\partial {x^{d]}}}{\partial {\alpha_n}}= \epsilon^{a b \cdots d} {\left\lvert{ \frac{\partial(x^1, x^2, \cdots x^n)}{\partial(\alpha_1, \alpha_2, \cdots \alpha_n)} }\right\rvert}\end{aligned} \hspace{\stretch{1}}(2.22)

When the dimension of the space $n$ is greater than the number of parameters for the integration hypervolume in question, the antisymmetric sum of partials is still the determinant of a Jacobian matrix

\begin{aligned}\frac{\partial {x^{[a_1}}}{\partial {\alpha_1}}\frac{\partial {x^{a_2}}}{\partial {\alpha_2}}\cdots \frac{\partial {x^{a_{n-1}}}}{\partial {\alpha_{n-1}}}\frac{\partial {x^{a_n]}}}{\partial {\alpha_n}}= {\left\lvert{ \frac{\partial(x^{a_1}, x^{a_2}, \cdots x^{a_n})}{\partial(\alpha_1, \alpha_2, \cdots \alpha_n)} }\right\rvert},\end{aligned} \hspace{\stretch{1}}(2.23)

however, we will have one such Jacobian for each unique choice of indexes.

# The Stokes work starts here.

The task is to relate our integral to the boundary of this volume, coming up with an explicit recipe for the description of that bounding surface, and determining the exact form of the reduced rank integral. This job is essentially to reduce the ranks of the tensors that are being contracted in our Stokes integral. With the derivative applied to our rank $n-1$ antisymmetric tensor $A_{b c \cdots d}$, we can apply the chain rule and examine the permutations so that this can be rewritten as a contraction of $A$ itself with a set of rank $n-1$ surface area elements.

\begin{aligned}\int d^n \alpha\frac{\partial {x^{[a}}}{\partial {\alpha_1}}\frac{\partial {x^b}}{\partial {\alpha_2}}\cdots \frac{\partial {x^c}}{\partial {\alpha_{n-1}}}\frac{\partial {x^{d]}}}{\partial {\alpha_n}}\partial_a A_{b c \cdots d} = ?\end{aligned} \hspace{\stretch{1}}(3.24)

Now, while the setup here has been completely general, this task is motivated by study of special relativity, where there is a requirement to work in a four dimensional space. Because of that explicit goal, I’m not going to attempt to formulate this in a completely abstract fashion. That task is really one of introducing sufficiently general notation. Instead, I’m going to proceed with a simpleton approach, and do this explicitly, and repeatedly for each of the rank 1, rank 2, and rank 3 tensor cases. It will be clear how this all generalizes by doing so, should one wish to work in still higher dimensional spaces.

## The rank 1 tensor case.

The equation we are working with for this vector case is

\begin{aligned}\int d^2 x \cdot (\nabla \wedge A) =\int d{\alpha_1} d{\alpha_2}\frac{\partial {x^{[a}}}{\partial {\alpha_1}}\frac{\partial {x^{b]}}}{\partial {\alpha_2}}\partial_a A_{b}(\alpha_1, \alpha_2)\end{aligned} \hspace{\stretch{1}}(3.25)

Expanding out the antisymmetric partials we have

\begin{aligned}\frac{\partial {x^{[a}}}{\partial {\alpha_1}}\frac{\partial {x^{b]}}}{\partial {\alpha_2}} & =\frac{\partial {x^{a}}}{\partial {\alpha_1}}\frac{\partial {x^{b}}}{\partial {\alpha_2}}-\frac{\partial {x^{b}}}{\partial {\alpha_1}}\frac{\partial {x^{a}}}{\partial {\alpha_2}},\end{aligned}

with which we can reduce the integral to

\begin{aligned}\int d^2 x \cdot (\nabla \wedge A) & =\int \left( d{\alpha_1}\frac{\partial {x^{a}}}{\partial {\alpha_1}}\frac{\partial {A_{b}}}{\partial {x^a}} \right)\frac{\partial {x^{b}}}{\partial {\alpha_2}} d{\alpha_2}-\left( d{\alpha_2}\frac{\partial {x^{a}}}{\partial {\alpha_2}}\frac{\partial {A_{b}}}{\partial {x^a}} \right)\frac{\partial {x^{b}}}{\partial {\alpha_1}} d{\alpha_1} \\ & =\int \left( d\alpha_1 \frac{\partial {A_b}}{\partial {\alpha_1}} \right)\frac{\partial {x^{b}}}{\partial {\alpha_2}} d{\alpha_2}-\left( d\alpha_2 \frac{\partial {A_b}}{\partial {\alpha_2}} \right)\frac{\partial {x^{b}}}{\partial {\alpha_1}} d{\alpha_1} \\ \end{aligned}

Now, if it happens that

\begin{aligned}\frac{\partial}{\partial {\alpha_1}}\frac{\partial {x^{a}}}{\partial {\alpha_2}} = \frac{\partial}{\partial {\alpha_2}}\frac{\partial {x^{a}}}{\partial {\alpha_1}} = 0\end{aligned} \hspace{\stretch{1}}(3.26)

then each of the individual integrals in $d\alpha_1$ and $d\alpha_2$ can be carried out. In that case, without any real loss of generality we can designate the integration bounds over the unit parametrization space square $\alpha_i \in [0,1]$, allowing this integral to be expressed as

\begin{aligned}\begin{aligned} & \int d{\alpha_1} d{\alpha_2}\frac{\partial {x^{[a}}}{\partial {\alpha_1}}\frac{\partial {x^{b]}}}{\partial {\alpha_2}}\partial_a A_{b}(\alpha_1, \alpha_2) \\ & =\int \left( A_b(1, \alpha_2) - A_b(0, \alpha_2) \right)\frac{\partial {x^{b}}}{\partial {\alpha_2}} d{\alpha_2}-\left( A_b(\alpha_1, 1) - A_b(\alpha_1, 0) \right)\frac{\partial {x^{b}}}{\partial {\alpha_1}} d{\alpha_1}.\end{aligned}\end{aligned} \hspace{\stretch{1}}(3.27)

It’s also fairly common to see ${\left.{{A}}\right\vert}_{{\partial \alpha_i}}$ used to designate evaluation of this first integral on the boundary, and using this we write

\begin{aligned}\int d{\alpha_1} d{\alpha_2}\frac{\partial {x^{[a}}}{\partial {\alpha_1}}\frac{\partial {x^{b]}}}{\partial {\alpha_2}}\partial_a A_{b}(\alpha_1, \alpha_2)=\int {\left.{{A_b}}\right\vert}_{{\partial \alpha_1}}\frac{\partial {x^{b}}}{\partial {\alpha_2}} d{\alpha_2}-{\left.{{A_b}}\right\vert}_{{\partial \alpha_2}}\frac{\partial {x^{b}}}{\partial {\alpha_1}} d{\alpha_1}.\end{aligned} \hspace{\stretch{1}}(3.28)

Also note that since we are summing over all $a,b$, and have

\begin{aligned}\frac{\partial {x^{[a}}}{\partial {\alpha_1}}\frac{\partial {x^{b]}}}{\partial {\alpha_2}}=-\frac{\partial {x^{[b}}}{\partial {\alpha_1}}\frac{\partial {x^{a]}}}{\partial {\alpha_2}},\end{aligned} \hspace{\stretch{1}}(3.29)

we can write this summing over all unique pairs of $a,b$ instead, which eliminates a small bit of redundancy (especially once the dimension of the vector space gets higher)

\begin{aligned}\boxed{\sum_{a < b}\int d{\alpha_1} d{\alpha_2}\frac{\partial {x^{[a}}}{\partial {\alpha_1}}\frac{\partial {x^{b]}}}{\partial {\alpha_2}}\left( \partial_a A_{b}-\partial_b A_{a} \right)=\int {\left.{{A_b}}\right\vert}_{{\partial \alpha_1}}\frac{\partial {x^{b}}}{\partial {\alpha_2}} d{\alpha_2}-{\left.{{A_b}}\right\vert}_{{\partial \alpha_2}}\frac{\partial {x^{b}}}{\partial {\alpha_1}} d{\alpha_1}.}\end{aligned} \hspace{\stretch{1}}(3.30)

In this form we have recovered the original geometric structure, with components of the curl multiplied by the component of the area element that shares the orientation and direction of that portion of the curl bivector.

This form of the result with evaluation at the boundaries in this form, assumed that ${\partial {x^a}}/{\partial {\alpha_1}}$ was not a function of $\alpha_2$ and ${\partial {x^a}}/{\partial {\alpha_2}}$ was not a function of $\alpha_1$. When that is not the case, we appear to have a less pretty result

\begin{aligned}\boxed{\sum_{a < b}\int d{\alpha_1} d{\alpha_2}\frac{\partial {x^{[a}}}{\partial {\alpha_1}}\frac{\partial {x^{b]}}}{\partial {\alpha_2}}\left( \partial_a A_{b}-\partial_b A_{a} \right)=\int d\alpha_2\int d\alpha_1\frac{\partial {A_b}}{\partial {\alpha_1}}\frac{\partial {x^{b}}}{\partial {\alpha_2}}-\int d\alpha_2\int d\alpha_1\frac{\partial {A_b}}{\partial {\alpha_2}}\frac{\partial {x^{b}}}{\partial {\alpha_1}}}\end{aligned} \hspace{\stretch{1}}(3.31)

Can this be reduced any further in the general case? Having seen the statements of Stokes theorem in it’s differential forms formulation, I initially expected the answer was yes, and only when I got to evaluating my $\mathbb{R}^{4}$ spacetime example below did I realize that the differentials displacements for the parallelogram that constituted the area element were functions of both parameters. Perhaps this detail is there in the differential forms version of the general Stokes theorem too, but is just hidden in a tricky fashion by the compact notation.

### Sanity check: $\mathbb{R}^{2}$ case in rectangular coordinates.

For $x^1 = x, x^2 = y$, and $\alpha_1 = x, \alpha_2 = y$, we have for the LHS

\begin{aligned} & \int_{x=x_0}^{x_1}\int_{y=y_0}^{y_1}dx dy\left(\frac{\partial {x^{1}}}{\partial {\alpha_1}}\frac{\partial {x^{2}}}{\partial {\alpha_2}}-\frac{\partial {x^{2}}}{\partial {\alpha_1}}\frac{\partial {x^{1}}}{\partial {\alpha_2}}\right)\partial_1 A_{2}+\left(\frac{\partial {x^{2}}}{\partial {\alpha_1}}\frac{\partial {x^{1}}}{\partial {\alpha_2}}-\frac{\partial {x^{1}}}{\partial {\alpha_1}}\frac{\partial {x^{2}}}{\partial {\alpha_2}}\right)\partial_2 A_{1} \\ & =\int_{x=x_0}^{x_1}\int_{y=y_0}^{y_1}dx dy\left( \frac{\partial {A_y}}{\partial x} - \frac{\partial {A_x}}{\partial y} \right)\end{aligned}

Our RHS expands to

\begin{aligned} & \int_{y=y_0}^{y_1} dy\left(\left( A_1(x_1, y) - A_1(x_0, y) \right)\frac{\partial {x^{1}}}{\partial y}+\left( A_2(x_1, y) - A_2(x_0, y) \right)\frac{\partial {x^{2}}}{\partial y}\right) \\ & \qquad-\int_{x=x_0}^{x_1} dx\left(\left( A_1(x, y_1) - A_1(x, y_0) \right)\frac{\partial {x^{1}}}{\partial x}+\left( A_2(x, y_1) - A_2(x, y_0) \right)\frac{\partial {x^{2}}}{\partial x}\right) \\ & =\int_{y=y_0}^{y_1} dy\left( A_y(x_1, y) - A_y(x_0, y) \right)-\int_{x=x_0}^{x_1} dx\left( A_x(x, y_1) - A_x(x, y_0) \right)\end{aligned}

We have

\begin{aligned}\begin{aligned} & \int_{x=x_0}^{x_1}\int_{y=y_0}^{y_1}dx dy\left( \frac{\partial {A_y}}{\partial x} - \frac{\partial {A_x}}{\partial y} \right) \\ & =\int_{y=y_0}^{y_1} dy\left( A_y(x_1, y) - A_y(x_0, y) \right)-\int_{x=x_0}^{x_1} dx\left( A_x(x, y_1) - A_x(x, y_0) \right)\end{aligned}\end{aligned} \hspace{\stretch{1}}(3.32)

The RHS is just a positively oriented line integral around the rectangle of integration

\begin{aligned}\int A_x(x, y_0) \hat{\mathbf{x}} \cdot ( \hat{\mathbf{x}} dx )+ A_y(x_1, y) \hat{\mathbf{y}} \cdot ( \hat{\mathbf{y}} dy )+ A_x(x, y_1) \hat{\mathbf{x}} \cdot ( -\hat{\mathbf{x}} dx )+ A_y(x_0, y) \hat{\mathbf{y}} \cdot ( -\hat{\mathbf{y}} dy )= \oint \mathbf{A} \cdot d\mathbf{r}.\end{aligned} \hspace{\stretch{1}}(3.33)

This special case is also recognizable as Green’s theorem, evident with the substitution $A_x = P$, $A_y = Q$, which gives us

\begin{aligned}\int_A dx dy \left( \frac{\partial Q}{\partial x} - \frac{\partial P}{\partial y} \right)=\oint_C P dx + Q dy.\end{aligned} \hspace{\stretch{1}}(3.34)

Strictly speaking, Green’s theorem is more general, since it applies to integration regions more general than rectangles, but that generalization can be arrived at easily enough, once the region is broken down into adjoining elementary regions.

### Sanity check: $\mathbb{R}^{3}$ case in rectangular coordinates.

It is expected that we can recover the classical Kelvin-Stokes theorem if we use rectangular coordinates in $\mathbb{R}^{3}$. However, we see that we have to consider three different parametrizations. If one picks rectangular parametrizations $(\alpha_1, \alpha_2) = \{ (x,y), (y,z), (z,x) \}$ in sequence, in each case holding the value of the additional coordinate fixed, we get three different independent Green’s function like relations

\begin{aligned}\int_A dx dy \left( \frac{\partial {A_y}}{\partial x} - \frac{\partial {A_x}}{\partial y} \right) & = \oint_C A_x dx + A_y dy \\ \int_A dy dz \left( \frac{\partial {A_z}}{\partial y} - \frac{\partial {A_y}}{\partial z} \right) & = \oint_C A_y dy + A_z dz \\ \int_A dz dx \left( \frac{\partial {A_x}}{\partial z} - \frac{\partial {A_z}}{\partial x} \right) & = \oint_C A_z dz + A_x dx.\end{aligned} \hspace{\stretch{1}}(3.35)

Note that we cannot just add these to form a complete integral $\oint \mathbf{A} \cdot d\mathbf{r}$ since the curves are all have different orientations. To recover the $\mathbb{R}^{3}$ Stokes theorem in rectangular coordinates, it appears that we’d have to consider a Riemann sum of triangular surface elements, and relate that to the loops over each of the surface elements. In that limiting argument, only the boundary of the complete surface would contribute to the RHS of the relation.

All that said, we shouldn’t actually have to go to all this work. Instead we can stick to a two variable parametrization of the surface, and use 3.30 directly.

### An illustration for a $\mathbb{R}^{4}$ spacetime surface.

Suppose we have a particle trajectory defined by an active Lorentz transformation from an initial spacetime point

\begin{aligned}x^i = O^{ij} x_j(0) = O^{ij} g_{jk} x^k = {O^{i}}_k x^k(0)\end{aligned} \hspace{\stretch{1}}(3.38)

Let the Lorentz transformation be formed by a composition of boost and rotation

\begin{aligned}{O^i}_j & = {L^i}_k {R^k}_j \\ {L^i}_j & =\begin{bmatrix}\cosh_\alpha & -\sinh\alpha & 0 & 0 \\ -\sinh_\alpha & \cosh\alpha & 0 & 0 \\ 0 & 0 & 1 & 0 \\ 0 & 0 & 0 & 1 \\ \end{bmatrix} \\ {R^i}_j & =\begin{bmatrix}1 & 0 & 0 & 0 \\ \cos_\alpha & \sin\alpha & 0 & 0 \\ -\sin_\alpha & \cos\alpha & 0 & 0 \\ 0 & 0 & 0 & 1 \\ \end{bmatrix}\end{aligned} \hspace{\stretch{1}}(3.39)

Different rates of evolution of $\alpha$ and $\theta$ define different trajectories, and taken together we have a surface described by the two parameters

\begin{aligned}x^i(\alpha, \theta) = {L^i}_k {R^k}_j x^j(0, 0).\end{aligned} \hspace{\stretch{1}}(3.42)

We can compute displacements along the trajectories formed by keeping either $\alpha$ or $\theta$ fixed and varying the other. Those are

\begin{aligned}\frac{\partial {x^i}}{\partial {\alpha}} d\alpha & = \frac{d{L^i}_k}{d\alpha} {R^k}_j x^j(0, 0) \\ \frac{\partial {x^i}}{\partial {\theta}} d\theta & = {L^i}_k \frac{d{R^k}_j}{d\theta} x^j(0, 0) .\end{aligned} \hspace{\stretch{1}}(3.43)

Writing $y^i = x^i(0,0)$ the computation of the partials above yields

\begin{aligned}\frac{\partial {x^i}}{\partial {\alpha}} & =\begin{bmatrix}\sinh\alpha y^0 -\cosh\alpha (\cos\theta y^1 + \sin\theta y^2) \\ -\cosh\alpha y^0 +\sinh\alpha (\cos\theta y^1 + \sin\theta y^2) \\ 0 \\ 0\end{bmatrix} \\ \frac{\partial {x^i}}{\partial {\theta}} & =\begin{bmatrix}-\sinh\alpha (-\sin\theta y^1 + \cos\theta y^2 ) \\ \cosh\alpha (-\sin\theta y^1 + \cos\theta y^2 ) \\ -(\cos\theta y^1 + \sin\theta y^2 ) \\ 0\end{bmatrix}\end{aligned} \hspace{\stretch{1}}(3.45)

Different choices of the initial point $y^i$ yield different surfaces, but we can get the idea by picking a simple starting point $y^i = (0, 1, 0, 0)$ leaving

\begin{aligned}\frac{\partial {x^i}}{\partial {\alpha}} & =\begin{bmatrix}-\cosh\alpha \cos\theta \\ \sinh\alpha \cos\theta \\ 0 \\ 0\end{bmatrix} \\ \frac{\partial {x^i}}{\partial {\theta}} & =\begin{bmatrix}\sinh\alpha \sin\theta \\ -\cosh\alpha \sin\theta \\ -\cos\theta \\ 0\end{bmatrix}.\end{aligned} \hspace{\stretch{1}}(3.47)

We can now compute our Jacobian determinants

\begin{aligned}\frac{\partial {x^{[a}}}{\partial {\alpha}} \frac{\partial {x^{b]}}}{\partial {\theta}}={\left\lvert{\frac{\partial(x^a, x^b)}{\partial(\alpha, \theta)}}\right\rvert}.\end{aligned} \hspace{\stretch{1}}(3.49)

Those are

\begin{aligned}{\left\lvert{\frac{\partial(x^0, x^1)}{\partial(\alpha, \theta)}}\right\rvert} & = \cos\theta \sin\theta \\ {\left\lvert{\frac{\partial(x^0, x^2)}{\partial(\alpha, \theta)}}\right\rvert} & = \cosh\alpha \cos^2\theta \\ {\left\lvert{\frac{\partial(x^0, x^3)}{\partial(\alpha, \theta)}}\right\rvert} & = 0 \\ {\left\lvert{\frac{\partial(x^1, x^2)}{\partial(\alpha, \theta)}}\right\rvert} & = -\sinh\alpha \cos^2\theta \\ {\left\lvert{\frac{\partial(x^1, x^3)}{\partial(\alpha, \theta)}}\right\rvert} & = 0 \\ {\left\lvert{\frac{\partial(x^2, x^3)}{\partial(\alpha, \theta)}}\right\rvert} & = 0\end{aligned} \hspace{\stretch{1}}(3.50)

Using this, let’s see a specific 4D example in spacetime for the integral of the curl of some four vector $A^i$, enumerating all the non-zero components of 3.31 for this particular spacetime surface

\begin{aligned}\sum_{a < b}\int d{\alpha} d{\theta}{\left\lvert{\frac{\partial(x^a, x^b)}{\partial(\alpha, \theta)}}\right\rvert}\left( \partial_a A_{b}-\partial_b A_{a} \right)=\int d\theta\int d\alpha\frac{\partial {A_b}}{\partial {\alpha}}\frac{\partial {x^{b}}}{\partial {\theta}}-\int d\theta\int d\alpha\frac{\partial {A_b}}{\partial {\theta}}\frac{\partial {x^{b}}}{\partial {\alpha}}\end{aligned} \hspace{\stretch{1}}(3.56)

The LHS is thus found to be

\begin{aligned} & \int d{\alpha} d{\theta}\left({\left\lvert{\frac{\partial(x^0, x^1)}{\partial(\alpha, \theta)}}\right\rvert} \left( \partial_0 A_{1} -\partial_1 A_{0} \right)+{\left\lvert{\frac{\partial(x^0, x^2)}{\partial(\alpha, \theta)}}\right\rvert} \left( \partial_0 A_{2} -\partial_2 A_{0} \right)+{\left\lvert{\frac{\partial(x^1, x^2)}{\partial(\alpha, \theta)}}\right\rvert} \left( \partial_1 A_{2} -\partial_2 A_{1} \right)\right) \\ & =\int d{\alpha} d{\theta}\left(\cos\theta \sin\theta \left( \partial_0 A_{1} -\partial_1 A_{0} \right)+\cosh\alpha \cos^2\theta \left( \partial_0 A_{2} -\partial_2 A_{0} \right)-\sinh\alpha \cos^2\theta \left( \partial_1 A_{2} -\partial_2 A_{1} \right)\right)\end{aligned}

On the RHS we have

\begin{aligned}\int d\theta\int d\alpha & \frac{\partial {A_b}}{\partial {\alpha}}\frac{\partial {x^{b}}}{\partial {\theta}}-\int d\theta\int d\alpha\frac{\partial {A_b}}{\partial {\theta}}\frac{\partial {x^{b}}}{\partial {\alpha}} \\ & =\int d\theta\int d\alpha\begin{bmatrix}\sinh\alpha \sin\theta & -\cosh\alpha \sin\theta & -\cos\theta & 0\end{bmatrix}\frac{\partial}{\partial {\alpha}}\begin{bmatrix}A_0 \\ A_1 \\ A_2 \\ A_3 \\ \end{bmatrix} \\ & -\int d\theta\int d\alpha\begin{bmatrix}-\cosh\alpha \cos\theta & \sinh\alpha \cos\theta & 0 & 0\end{bmatrix}\frac{\partial}{\partial {\theta}}\begin{bmatrix}A_0 \\ A_1 \\ A_2 \\ A_3 \\ \end{bmatrix} \\ \end{aligned}

\begin{aligned}\begin{aligned} & \int d{\alpha} d{\theta}\cos\theta \sin\theta \left( \partial_0 A_{1} -\partial_1 A_{0} \right) \\ & \qquad+\int d{\alpha} d{\theta}\cosh\alpha \cos^2\theta \left( \partial_0 A_{2} -\partial_2 A_{0} \right) \\ & \qquad-\int d{\alpha} d{\theta}\sinh\alpha \cos^2\theta \left( \partial_1 A_{2} -\partial_2 A_{1} \right) \\ & =\int d\theta \sin\theta \int d\alpha \left( \sinh\alpha \frac{\partial {A_0}}{\partial {\alpha}} - \cosh\alpha \frac{\partial {A_1}}{\partial {\alpha}} \right) \\ & \qquad-\int d\theta \cos\theta \int d\alpha \frac{\partial {A_2}}{\partial {\alpha}} \\ & \qquad+\int d\alpha \cosh\alpha \int d\theta \cos\theta \frac{\partial {A_0}}{\partial {\theta}} \\ & \qquad-\int d\alpha \sinh\alpha \int d\theta \cos\theta \frac{\partial {A_1}}{\partial {\theta}}\end{aligned}\end{aligned} \hspace{\stretch{1}}(3.57)

Because of the complexity of the surface, only the second term on the RHS has the “evaluate on the boundary” characteristic that may have been expected from a Green’s theorem like line integral.

It is also worthwhile to point out that we have had to be very careful with upper and lower indexes all along (and have done so with the expectation that our application would include the special relativity case where our metric determinant is minus one.) Because we worked with upper indexes for the area element, we had to work with lower indexes for the four vector and the components of the gradient that we included in our curl evaluation.

## The rank 2 tensor case.

Let’s consider briefly the terms in the contraction sum

\begin{aligned}{\left\lvert{ \frac{\partial(x^a, x^b, x^c)}{\partial(\alpha_1, \alpha_2, \alpha_3)} }\right\rvert} \partial_a A_{bc}\end{aligned} \hspace{\stretch{1}}(3.58)

For any choice of a set of three distinct indexes $(a, b, c) \in (0, 1, 2), (0, 1, 3), (0, 2, 3), (1, 2, 3)$), we have $6 = 3!$ ways of permuting those indexes in this sum

\begin{aligned}{\left\lvert{ \frac{\partial(x^a, x^b, x^c)}{\partial(\alpha_1, \alpha_2, \alpha_3)} }\right\rvert} \partial_a A_{bc} & =\sum_{a < b < c} {\left\lvert{ \frac{\partial(x^a, x^b, x^c)}{\partial(\alpha_1, \alpha_2, \alpha_3)} }\right\rvert} \partial_a A_{bc} + {\left\lvert{ \frac{\partial(x^a, x^c, x^b)}{\partial(\alpha_1, \alpha_2, \alpha_3)} }\right\rvert} \partial_a A_{cb} + {\left\lvert{ \frac{\partial(x^b, x^c, x^a)}{\partial(\alpha_1, \alpha_2, \alpha_3)} }\right\rvert} \partial_b A_{ca} \\ & \qquad + {\left\lvert{ \frac{\partial(x^b, x^a, x^c)}{\partial(\alpha_1, \alpha_2, \alpha_3)} }\right\rvert} \partial_b A_{ac} + {\left\lvert{ \frac{\partial(x^c, x^a, x^b)}{\partial(\alpha_1, \alpha_2, \alpha_3)} }\right\rvert} \partial_c A_{ab} + {\left\lvert{ \frac{\partial(x^c, x^b, x^a)}{\partial(\alpha_1, \alpha_2, \alpha_3)} }\right\rvert} \partial_c A_{ba} \\ & =2!\sum_{a < b < c}{\left\lvert{ \frac{\partial(x^a, x^b, x^c)}{\partial(\alpha_1, \alpha_2, \alpha_3)} }\right\rvert}\left( \partial_a A_{bc} + \partial_b A_{c a} + \partial_c A_{a b} \right)\end{aligned}

Observe that we have no sign alternation like we had in the vector (rank 1 tensor) case. That sign alternation in this summation expansion appears to occur only for odd grade tensors.

Returning to the problem, we wish to expand the determinant in order to apply a chain rule contraction as done in the rank-1 case. This can be done along any of rows or columns of the determinant, and we can write any of

\begin{aligned}{\left\lvert{ \frac{\partial(x^a, x^b, x^c)}{\partial(\alpha_1, \alpha_2, \alpha_3)} }\right\rvert} & =\frac{\partial {x^a}}{\partial {\alpha_1}} {\left\lvert{ \frac{\partial(x^b, x^c)}{\partial(\alpha_2, \alpha_3)} }\right\rvert}-\frac{\partial {x^a}}{\partial {\alpha_2}} {\left\lvert{ \frac{\partial(x^b, x^c)}{\partial(\alpha_1, \alpha_3)} }\right\rvert}+\frac{\partial {x^a}}{\partial {\alpha_3}} {\left\lvert{ \frac{\partial(x^b, x^c)}{\partial(\alpha_1, \alpha_2)} }\right\rvert} \\ & =\frac{\partial {x^b}}{\partial {\alpha_1}} {\left\lvert{ \frac{\partial(x^c, x^a)}{\partial(\alpha_2, \alpha_3)} }\right\rvert}-\frac{\partial {x^b}}{\partial {\alpha_2}} {\left\lvert{ \frac{\partial(x^c, x^a)}{\partial(\alpha_1, \alpha_3)} }\right\rvert}+\frac{\partial {x^b}}{\partial {\alpha_3}} {\left\lvert{ \frac{\partial(x^c, x^a)}{\partial(\alpha_1, \alpha_2)} }\right\rvert} \\ & =\frac{\partial {x^c}}{\partial {\alpha_1}} {\left\lvert{ \frac{\partial(x^a, x^b)}{\partial(\alpha_2, \alpha_3)} }\right\rvert}-\frac{\partial {x^c}}{\partial {\alpha_2}} {\left\lvert{ \frac{\partial(x^a, x^b)}{\partial(\alpha_1, \alpha_3)} }\right\rvert}+\frac{\partial {x^c}}{\partial {\alpha_3}} {\left\lvert{ \frac{\partial(x^a, x^b)}{\partial(\alpha_1, \alpha_2)} }\right\rvert} \\ \end{aligned}

This allows the contraction of the index $a$, eliminating it from the result

\begin{aligned}{\left\lvert{ \frac{\partial(x^a, x^b, x^c)}{\partial(\alpha_1, \alpha_2, \alpha_3)} }\right\rvert} \partial_a A_{bc} & =\left( \frac{\partial {x^a}}{\partial {\alpha_1}} {\left\lvert{ \frac{\partial(x^b, x^c)}{\partial(\alpha_2, \alpha_3)} }\right\rvert}-\frac{\partial {x^a}}{\partial {\alpha_2}} {\left\lvert{ \frac{\partial(x^b, x^c)}{\partial(\alpha_1, \alpha_3)} }\right\rvert}+\frac{\partial {x^a}}{\partial {\alpha_3}} {\left\lvert{ \frac{\partial(x^b, x^c)}{\partial(\alpha_1, \alpha_2)} }\right\rvert} \right) \frac{\partial {A_{bc}}}{\partial {x^a}} \\ & =\frac{\partial {A_{bc}}}{\partial {\alpha_1}} {\left\lvert{ \frac{\partial(x^b, x^c)}{\partial(\alpha_2, \alpha_3)} }\right\rvert}-\frac{\partial {A_{bc}}}{\partial {\alpha_2}} {\left\lvert{ \frac{\partial(x^b, x^c)}{\partial(\alpha_1, \alpha_3)} }\right\rvert}+\frac{\partial {A_{bc}}}{\partial {\alpha_3}} {\left\lvert{ \frac{\partial(x^b, x^c)}{\partial(\alpha_1, \alpha_2)} }\right\rvert} \\ & =2!\sum_{b < c}\frac{\partial {A_{bc}}}{\partial {\alpha_1}} {\left\lvert{ \frac{\partial(x^b, x^c)}{\partial(\alpha_2, \alpha_3)} }\right\rvert}-\frac{\partial {A_{bc}}}{\partial {\alpha_2}} {\left\lvert{ \frac{\partial(x^b, x^c)}{\partial(\alpha_1, \alpha_3)} }\right\rvert}+\frac{\partial {A_{bc}}}{\partial {\alpha_3}} {\left\lvert{ \frac{\partial(x^b, x^c)}{\partial(\alpha_1, \alpha_2)} }\right\rvert} \\ \end{aligned}

Dividing out the common $2!$ terms, we can summarize this result as

\begin{aligned}\boxed{\begin{aligned}\sum_{a < b < c} & \int d\alpha_1 d\alpha_2 d\alpha_3 {\left\lvert{ \frac{\partial(x^a, x^b, x^c)}{\partial(\alpha_1, \alpha_2, \alpha_3)} }\right\rvert}\left( \partial_a A_{bc} + \partial_b A_{c a} + \partial_c A_{a b} \right) \\ & =\sum_{b < c}\int d\alpha_2 d\alpha_3 \int d\alpha_1\frac{\partial {A_{bc}}}{\partial {\alpha_1}} {\left\lvert{ \frac{\partial(x^b, x^c)}{\partial(\alpha_2, \alpha_3)} }\right\rvert} \\ & -\sum_{b < c}\int d\alpha_1 d\alpha_3 \int d\alpha_2\frac{\partial {A_{bc}}}{\partial {\alpha_2}} {\left\lvert{ \frac{\partial(x^b, x^c)}{\partial(\alpha_1, \alpha_3)} }\right\rvert} \\ & +\sum_{b < c}\int d\alpha_1 d\alpha_2 \int d\alpha_3\frac{\partial {A_{bc}}}{\partial {\alpha_3}} {\left\lvert{ \frac{\partial(x^b, x^c)}{\partial(\alpha_1, \alpha_2)} }\right\rvert}\end{aligned}}\end{aligned} \hspace{\stretch{1}}(3.59)

In general, as observed in the spacetime surface example above, the two index Jacobians can be functions of the integration variable first being eliminated. In the special cases where this is not the case (such as the $\mathbb{R}^{3}$ case with rectangular coordinates), then we are left with just the evaluation of the tensor element $A_{bc}$ on the boundaries of the respective integrals.

## The rank 3 tensor case.

The key step is once again just a determinant expansion

\begin{aligned} {\partial(\alpha_1, \alpha_2, \alpha_3, \alpha_4)} }\right\rvert} \\ & =\frac{\partial {x^a}}{\partial {\alpha_1}} {\left\lvert{ \frac{\partial(x^b, x^c, x^d)}{\partial(\alpha_2, \alpha_3, \alpha_4)} }\right\rvert}-\frac{\partial {x^a}}{\partial {\alpha_2}} {\left\lvert{ \frac{\partial(x^b, x^c, x^d)}{\partial(\alpha_1, \alpha_3, \alpha_4)} }\right\rvert}+\frac{\partial {x^a}}{\partial {\alpha_3}} {\left\lvert{ \frac{\partial(x^b, x^c, x^d)}{\partial(\alpha_1, \alpha_2, \alpha_4)} }\right\rvert}+\frac{\partial {x^a}}{\partial {\alpha_4}} {\left\lvert{ \frac{\partial(x^b, x^c, x^d)}{\partial(\alpha_1, \alpha_2, \alpha_3)} }\right\rvert}\\ \end{aligned}

so that the sum can be reduced from a four index contraction to a 3 index contraction

\begin{aligned} {\partial(\alpha_1, \alpha_2, \alpha_3, \alpha_4)} }\right\rvert} \partial_a A_{bcd} \\ & =\frac{\partial {A_{bcd}}}{\partial {\alpha_1}} {\left\lvert{ \frac{\partial(x^b, x^c, x^d)}{\partial(\alpha_2, \alpha_3, \alpha_4)} }\right\rvert}-\frac{\partial {A_{bcd}}}{\partial {\alpha_2}} {\left\lvert{ \frac{\partial(x^b, x^c, x^d)}{\partial(\alpha_1, \alpha_3, \alpha_4)} }\right\rvert}+\frac{\partial {A_{bcd}}}{\partial {\alpha_3}} {\left\lvert{ \frac{\partial(x^b, x^c, x^d)}{\partial(\alpha_1, \alpha_2, \alpha_4)} }\right\rvert}+\frac{\partial {A_{bcd}}}{\partial {\alpha_4}} {\left\lvert{ \frac{\partial(x^b, x^c, x^d)}{\partial(\alpha_1, \alpha_2, \alpha_3)} }\right\rvert}\end{aligned}

That’s the essence of the theorem, but we can play the same combinatorial reduction games to reduce the built in redundancy in the result

\begin{aligned}\boxed{\begin{aligned}\frac{1}{{3!}} & \int d^4 \alpha {\left\lvert{ \frac{\partial(x^a, x^b, x^c, x^d)}{\partial(\alpha_1, \alpha_2, \alpha_3, \alpha_4)} }\right\rvert} \partial_a A_{bcd} \\ & =\sum_{a < b < c < d}\int d^4 \alpha {\left\lvert{ \frac{\partial(x^a, x^b, x^c, x^d)}{\partial(\alpha_1, \alpha_2, \alpha_3, \alpha_4)} }\right\rvert} \left( \partial_a A_{bcd} -\partial_b A_{cda} +\partial_c A_{dab} -\partial_d A_{abc} \right) \\ & =\qquad \sum_{b < c < d}\int d\alpha_2 d\alpha_3 d\alpha_4 \int d\alpha_1\frac{\partial {A_{bcd}}}{\partial {\alpha_1}} {\left\lvert{ \frac{\partial(x^b, x^c, x^d)}{\partial(\alpha_2, \alpha_3, \alpha_4)} }\right\rvert} \\ & \qquad -\sum_{b < c < d}\int d\alpha_1 d\alpha_3 d\alpha_4 \int d\alpha_2\frac{\partial {A_{bcd}}}{\partial {\alpha_2}} {\left\lvert{ \frac{\partial(x^b, x^c, x^d)}{\partial(\alpha_1, \alpha_3, \alpha_4)} }\right\rvert} \\ & \qquad +\sum_{b < c < d}\int d\alpha_1 d\alpha_2 d\alpha_4 \int d\alpha_3\frac{\partial {A_{bcd}}}{\partial {\alpha_3}} {\left\lvert{ \frac{\partial(x^b, x^c, x^d)}{\partial(\alpha_1, \alpha_2, \alpha_4)} }\right\rvert} \\ & \qquad +\sum_{b < c < d}\int d\alpha_1 d\alpha_2 d\alpha_3 \int d\alpha_4\frac{\partial {A_{bcd}}}{\partial {\alpha_4}} {\left\lvert{ \frac{\partial(x^b, x^c, x^d)}{\partial(\alpha_1, \alpha_2, \alpha_3)} }\right\rvert} \\ \end{aligned}}\end{aligned} \hspace{\stretch{1}}(3.60)

## A note on Four diverence.

Our four divergence integral has the following form

\begin{aligned}\int d^4 \alpha {\left\lvert{ \frac{\partial(x^1, x^2, x^2, x^4)}{\partial(\alpha_1, \alpha_2, \alpha_3, \alpha_4)} }\right\rvert} \partial_a A^a\end{aligned} \hspace{\stretch{1}}(3.61)

We can relate this to the rank 3 Stokes theorem with a duality transformation, multiplying with a pseudoscalar

\begin{aligned}A^a = \epsilon^{abcd} T_{bcd},\end{aligned} \hspace{\stretch{1}}(3.62)

where $T_{bcd}$ can also be related back to the vector by the same sort of duality transformation

\begin{aligned}A^a \epsilon_{a b c d} = \epsilon^{abcd} \epsilon_{a b c d} T_{bcd} = 4! T_{bcd}.\end{aligned} \hspace{\stretch{1}}(3.63)

The divergence integral in terms of the rank 3 tensor is

\begin{aligned}\int d^4 \alpha {\left\lvert{ \frac{\partial(x^1, x^2, x^2, x^4)}{\partial(\alpha_1, \alpha_2, \alpha_3, \alpha_4)} }\right\rvert} \partial_a \epsilon^{abcd} T_{bcd}=\int d^4 \alpha {\left\lvert{ \frac{\partial(x^a, x^b, x^c, x^d)}{\partial(\alpha_1, \alpha_2, \alpha_3, \alpha_4)} }\right\rvert} \partial_a T_{bcd},\end{aligned} \hspace{\stretch{1}}(3.64)

and we are free to perform the same Stokes reduction of the integral. Of course, this is particularly simple in rectangular coordinates. I still have to think though one sublty that I feel may be important. We could have started off with an integral of the following form

\begin{aligned}\int dx^1 dx^2 dx^3 dx^4 \partial_a A^a,\end{aligned} \hspace{\stretch{1}}(3.65)

and I think this differs from our starting point slightly because this has none of the antisymmetric structure of the signed 4 volume element that we have used. We do not take the absolute value of our Jacobians anywhere.

## Jacobians and spherical polar gradient

Posted by peeterjoot on December 8, 2009

[Click here for a PDF of this post with nicer formatting]

# Motivation

The dumbest and most obvious way to do a chain of variables for the gradient is to utilize a chain rule expansion producing the Jacobian matrix to transform the coordinates. Here we do this to calculate the spherical polar representation of the gradient.

There are smarter and easier ways to do this, but there is some surprising simple structure to the resulting Jacobians that seems worth noting.

# Spherical polar gradient coordinates in terms of Cartesian.

We wish to do a change of variables for each of the differential operators of the gradient. This is essentially just application of the chain rule, as in

\begin{aligned}\frac{\partial {}}{\partial {r}} = \frac{\partial {x}}{\partial {r}} \frac{\partial {}}{\partial {x}}+\frac{\partial {y}}{\partial {r}} \frac{\partial {}}{\partial {y}}+\frac{\partial {z}}{\partial {r}} \frac{\partial {}}{\partial {z}}.\end{aligned} \quad\quad\quad(1)

Collecting all such derivatives we have in column vector form

\begin{aligned}\begin{bmatrix}\partial_r \\ \partial_\theta \\ \partial_\phi\end{bmatrix}= \begin{bmatrix}\frac{\partial {x}}{\partial {r}} &\frac{\partial {y}}{\partial {r}} &\frac{\partial {z}}{\partial {r}} \\ \frac{\partial {x}}{\partial {\theta}} &\frac{\partial {y}}{\partial {\theta}} &\frac{\partial {z}}{\partial {\theta}} \\ \frac{\partial {x}}{\partial {\phi}} &\frac{\partial {y}}{\partial {\phi}} &\frac{\partial {z}}{\partial {\phi}} \end{bmatrix}\begin{bmatrix}\partial_x \\ \partial_y \\ \partial_z\end{bmatrix}.\end{aligned} \quad\quad\quad(2)

This becomes a bit more tractable with the Jacobian notation

\begin{aligned}\frac{\partial (x,y,z)}{\partial (r,\theta,\phi)}=\begin{bmatrix}\frac{\partial {x}}{\partial {r}} &\frac{\partial {y}}{\partial {r}} &\frac{\partial {z}}{\partial {r}} \\ \frac{\partial {x}}{\partial {\theta}} &\frac{\partial {y}}{\partial {\theta}} &\frac{\partial {z}}{\partial {\theta}} \\ \frac{\partial {x}}{\partial {\phi}} &\frac{\partial {y}}{\partial {\phi}} &\frac{\partial {z}}{\partial {\phi}}\end{bmatrix}.\end{aligned} \quad\quad\quad(3)

The change of variables for the operator triplet is then just

\begin{aligned}\begin{bmatrix}\partial_r \\ \partial_\theta \\ \partial_\phi\end{bmatrix}= \frac{\partial (x,y,z)}{\partial (r,\theta,\phi)}\begin{bmatrix}\partial_x \\ \partial_y \\ \partial_z\end{bmatrix}.\end{aligned} \quad\quad\quad(4)

This Jacobian matrix is also not even too hard to calculate. With $\mathbf{x} = r \hat{\mathbf{r}}$, we have $x_k = r \hat{\mathbf{r}} \cdot \mathbf{e}_k$, and

\begin{aligned}\frac{\partial {x_k}}{\partial {r}} &= \hat{\mathbf{r}} \cdot \mathbf{e}_k \\ \frac{\partial {x_k}}{\partial {\theta}} &= r \frac{\partial {\hat{\mathbf{r}}}}{\partial {\theta}} \cdot \mathbf{e}_k \\ \frac{\partial {x_k}}{\partial {\phi}} &= r \frac{\partial {\hat{\mathbf{r}}}}{\partial {\phi}} \cdot \mathbf{e}_k.\end{aligned} \quad\quad\quad(5)

The last two derivatives can be calculated easily if the radial unit vector is written out explicitly, with $S$ and $C$ for sine and cosine respectively, these are

\begin{aligned}\hat{\mathbf{r}} &= \begin{bmatrix}S_\theta C_\phi \\ S_\theta S_\phi \\ C_\theta \end{bmatrix} \\ \frac{\partial {\hat{\mathbf{r}}}}{\partial {\theta}} &= \begin{bmatrix}C_\theta C_\phi \\ C_\theta S_\phi \\ -S_\theta \end{bmatrix} \\ \frac{\partial {\hat{\mathbf{r}}}}{\partial {\phi}} &= \begin{bmatrix}-S_\theta S_\phi \\ S_\theta C_\phi \\ 0\end{bmatrix} .\end{aligned} \quad\quad\quad(8)

We can plug these into the elements of the Jacobian matrix explicitly, which produces

\begin{aligned}\frac{\partial (x,y,z)}{\partial (r,\theta,\phi)}=\begin{bmatrix} S_\theta C_\phi & S_\theta S_\phi & C_\theta \\ r C_\theta C_\phi & r C_\theta S_\phi & - r S_\theta \\ -r S_\theta S_\phi & rS_\theta C_\phi & 0\end{bmatrix},\end{aligned} \quad\quad\quad(11)

however, we are probably better off just referring back to 8, and writing

\begin{aligned}\frac{\partial (x,y,z)}{\partial (r,\theta,\phi)}=\begin{bmatrix} \hat{\mathbf{r}}^\text{T} \\ r \frac{\partial {\hat{\mathbf{r}}^\text{T}}}{\partial {\theta}} \\ r \frac{\partial {\hat{\mathbf{r}}^\text{T}}}{\partial {\phi}} \end{bmatrix}.\end{aligned} \quad\quad\quad(12)

Unfortunately, this is actually a bit of a dead end. We really want the inverse of this matrix because the desired quantity is

\begin{aligned}\boldsymbol{\nabla} = \begin{bmatrix}\mathbf{e}_1 & \mathbf{e}_2 & \mathbf{e}_3 \end{bmatrix}\begin{bmatrix}\partial_{x_1} \\ \partial_{x_2} \\ \partial_{x_3}\end{bmatrix}.\end{aligned} \quad\quad\quad(13)

(Here my matrix of unit vectors treats these abusively as single elements and not as column vectors).

The matrix of equation 12 does not look particularly fun to invert directly, and that is what we need to substitute into
13. One knows that in the end if it was attempted things should mystically simplify (presuming this was done error free).

# Cartesian gradient coordinates in terms of spherical polar partials.

Let’s flip things upside down and calculate the inverse Jacobian matrix directly. This is a messier job, but it appears less messy than the matrix inversion above.

\begin{aligned}r^2 &= x^2 + y^2 + z^2 \\ \sin^2 \theta &= \frac{x^2 + y^2}{x^2 + y^2 + z^2} \\ \tan\phi &= \frac{y}{x}.\end{aligned} \quad\quad\quad(14)

The messy task is now the calculation of these derivatives.

For the first, from $r^2 = x^2 + y^2 + z^2$, taking partials on both sides, we have

\begin{aligned}\frac{\partial {r}}{\partial {x_k}} = \frac{x_k}{r}.\end{aligned} \quad\quad\quad(17)

But these are just the direction cosines, the components of our polar unit vector $\hat{\mathbf{r}}$. We can then write for all of these derivatives in column matrix form

\begin{aligned}\boldsymbol{\nabla} r = \hat{\mathbf{r}}\end{aligned} \quad\quad\quad(18)

Next from $\sin^2\theta = (x^2 + y^2)/r^2$, we get after some reduction

\begin{aligned}\frac{\partial {\theta}}{\partial {x}} &= \frac{1}{{r}} C_\theta C_\phi \\ \frac{\partial {\phi}}{\partial {y}} &= \frac{1}{{r}} C_\theta S_\phi \\ \frac{\partial {\phi}}{\partial {z}} &= -\frac{S_\theta}{r}.\end{aligned} \quad\quad\quad(19)

Observe that we can antidifferentiate with respect to theta and obtain

\begin{aligned}\boldsymbol{\nabla} \theta &= \frac{1}{{r}}\begin{bmatrix}C_\theta C_\phi \\ C_\theta S_\phi \\ -S_\theta\end{bmatrix} \\ &=\frac{1}{{r}}\frac{\partial {}}{\partial {\theta}}\begin{bmatrix}S_\theta C_\phi \\ S_\theta S_\phi \\ C_\theta\end{bmatrix}.\end{aligned}

This last column vector is our friend the unit polar vector again, and we have

\begin{aligned}\boldsymbol{\nabla} \theta &= \frac{1}{{r}}\frac{\partial {\hat{\mathbf{r}}}}{\partial {\theta}}\end{aligned} \quad\quad\quad(22)

Finally for the $\phi$ dependence we have after some reduction

\begin{aligned}\boldsymbol{\nabla} \phi &=\frac{1}{{r S_\theta}}\begin{bmatrix}-S_\phi \\ C_\phi \\ 0\end{bmatrix}.\end{aligned} \quad\quad\quad(23)

Again, we can antidifferentiate

\begin{aligned}\boldsymbol{\nabla} \phi &=\frac{1}{{r (S_\theta)^2}}\begin{bmatrix}-S_\theta S_\phi \\ S_\theta C_\phi \\ 0\end{bmatrix} \\ &=\frac{1}{{r (S_\theta)^2}}\frac{\partial {}}{\partial {\phi}}\begin{bmatrix}S_\theta C_\phi \\ S_\theta S_\phi \\ C_\theta\end{bmatrix}.\end{aligned}

We have our unit polar vector again, and our $\phi$ partials nicely summarized by

\begin{aligned}\boldsymbol{\nabla} \phi &=\frac{1}{{r (S_\theta)^2}}\frac{\partial {\hat{\mathbf{r}}}}{\partial {\phi}}.\end{aligned} \quad\quad\quad(24)

With this we can now write out the Jacobian matrix either explicitly, or in column vector form in terms of $\hat{\mathbf{r}}$. First a reminder of why we want this matrix, for the following change of variables

\begin{aligned}\begin{bmatrix}\partial_x \\ \partial_y \\ \partial_z\end{bmatrix}= \begin{bmatrix}\frac{\partial {r}}{\partial {x}} &\frac{\partial {\theta}}{\partial {x}} &\frac{\partial {\phi}}{\partial {x}} \\ \frac{\partial {r}}{\partial {y}} &\frac{\partial {\theta}}{\partial {y}} &\frac{\partial {\phi}}{\partial {y}} \\ \frac{\partial {r}}{\partial {z}} &\frac{\partial {\theta}}{\partial {z}} &\frac{\partial {\phi}}{\partial {z}} \end{bmatrix}\begin{bmatrix}\partial_r \\ \partial_\theta \\ \partial_\phi\end{bmatrix}.\end{aligned} \quad\quad\quad(25)

We want the Jacobian matrix

\begin{aligned}\frac{\partial (r,\theta,\phi)}{\partial (x, y, z)}=\begin{bmatrix}\boldsymbol{\nabla} r & \boldsymbol{\nabla} \theta & \boldsymbol{\nabla} \phi\end{bmatrix}=\begin{bmatrix}\hat{\mathbf{r}} & \frac{1}{{r}} \frac{\partial {\hat{\mathbf{r}}}}{\partial {\theta}} & \frac{1}{{r \sin^2\theta}} \frac{\partial {\hat{\mathbf{r}}}}{\partial {\phi}}\end{bmatrix}.\end{aligned} \quad\quad\quad(26)

Explicitly this is

\begin{aligned}\frac{\partial (r,\theta,\phi)}{\partial (x, y, z)}=\begin{bmatrix}S_\theta C_\phi & \frac{1}{{r}} C_\theta C_\phi & -\frac{1}{{r S_\theta}} S_\phi \\ S_\theta S_\phi & \frac{1}{{r}} C_\theta S_\phi & \frac{C_\phi}{r S_\theta} \\ C_\theta & -\frac{1}{{r}} S_\theta & 0\end{bmatrix}.\end{aligned} \quad\quad\quad(27)

As a verification of correctness multiplication of this with 11 should produce identity. That’s a mess of trig that I don’t really feel like trying, but we can get a rough idea why it should all be the identity matrix by multiplying it out in block matrix form

\begin{aligned}\frac{\partial (x,y,z)}{\partial (r,\theta,\phi)}\frac{\partial (r,\theta,\phi)}{\partial (x, y, z)}&=\begin{bmatrix} \hat{\mathbf{r}}^\text{T} \\ r \frac{\partial {\hat{\mathbf{r}}^\text{T}}}{\partial {\theta}} \\ r \frac{\partial {\hat{\mathbf{r}}^\text{T}}}{\partial {\phi}} \end{bmatrix}\begin{bmatrix}\hat{\mathbf{r}} & \frac{1}{{r}} \frac{\partial {\hat{\mathbf{r}}}}{\partial {\theta}} & \frac{1}{{r \sin^2\theta}} \frac{\partial {\hat{\mathbf{r}}}}{\partial {\phi}}\end{bmatrix} \\ &=\begin{bmatrix} \hat{\mathbf{r}}^\text{T} \hat{\mathbf{r}} & \frac{1}{{r}} \hat{\mathbf{r}}^\text{T} \frac{\partial {\hat{\mathbf{r}}}}{\partial {\theta}} & \frac{1}{{r \sin^2 \theta}} \hat{\mathbf{r}}^\text{T} \frac{\partial {\hat{\mathbf{r}}}}{\partial {\phi}} \\ r \frac{\partial {\hat{\mathbf{r}}^\text{T}}}{\partial {\theta}} \hat{\mathbf{r}} & \frac{\partial {\hat{\mathbf{r}}^\text{T}}}{\partial {\theta}} \frac{\partial {\hat{\mathbf{r}}}}{\partial {\theta}} & \frac{1}{{\sin^2\theta}} \frac{\partial {\hat{\mathbf{r}}^\text{T}}}{\partial {\theta}} \frac{\partial {\hat{\mathbf{r}}}}{\partial {\phi}} \\ r \frac{\partial {\hat{\mathbf{r}}^\text{T}}}{\partial {\phi}} \hat{\mathbf{r}} & \frac{\partial {\hat{\mathbf{r}}^\text{T}}}{\partial {\phi}} \frac{\partial {\hat{\mathbf{r}}}}{\partial {\theta}} & \frac{1}{{\sin^2\theta}} \frac{\partial {\hat{\mathbf{r}}^\text{T}}}{\partial {\phi}} \frac{\partial {\hat{\mathbf{r}}}}{\partial {\phi}}\end{bmatrix}.\end{aligned}

The derivatives are vectors that lie tangential to the unit sphere. We can calculate this to verify, or we can look at the off diagonal terms which say just this if we trust the math that says these should all be zeros. For each of the off diagonal terms to be zero must mean that we have

\begin{aligned}0 = \frac{\partial {\hat{\mathbf{r}}}}{\partial {\theta}} \cdot \hat{\mathbf{r}} = \frac{\partial {\hat{\mathbf{r}}}}{\partial {\theta}} \cdot \frac{\partial {\hat{\mathbf{r}}}}{\partial {\phi}} = \frac{\partial {\hat{\mathbf{r}}}}{\partial {\phi}} \cdot \hat{\mathbf{r}} \end{aligned} \quad\quad\quad(28)

This makes intuitive sense. We can also verify quickly enough that $({\partial {\hat{\mathbf{r}}}}/{\partial {\theta}})^2 = 1$, and $({\partial {\hat{\mathbf{r}}}}/{\partial {\phi}})^2 = \sin^2\theta$ (I did this with a back of the envelope calculation using geometric algebra). That is consistent with what this matrix product implies it should equal.

# Completing the gradient change of variables to spherical polar coordinates.

We are now set to calculate the gradient in spherical polar coordinates from our Cartesian representation. From 13 and
25, and 26 we have

\begin{aligned}\boldsymbol{\nabla} =\begin{bmatrix}\mathbf{e}_1 & \mathbf{e}_2 & \mathbf{e}_3 \end{bmatrix}\begin{bmatrix}\hat{\mathbf{r}} \cdot \mathbf{e}_1 & \frac{1}{{r}} \frac{\partial {\hat{\mathbf{r}}}}{\partial {\theta}} \cdot \mathbf{e}_1 & \frac{1}{{r \sin^2\theta}} \frac{\partial {\hat{\mathbf{r}}}}{\partial {\phi}} \cdot \mathbf{e}_1 \\ \hat{\mathbf{r}} \cdot \mathbf{e}_2 & \frac{1}{{r}} \frac{\partial {\hat{\mathbf{r}}}}{\partial {\theta}} \cdot \mathbf{e}_2 & \frac{1}{{r \sin^2\theta}} \frac{\partial {\hat{\mathbf{r}}}}{\partial {\phi}} \cdot \mathbf{e}_2 \\ \hat{\mathbf{r}} \cdot \mathbf{e}_3 & \frac{1}{{r}} \frac{\partial {\hat{\mathbf{r}}}}{\partial {\theta}} \cdot \mathbf{e}_3 & \frac{1}{{r \sin^2\theta}} \frac{\partial {\hat{\mathbf{r}}}}{\partial {\phi}} \cdot \mathbf{e}_3 \end{bmatrix}\begin{bmatrix}\partial_r \\ \partial_\theta \\ \partial_\phi\end{bmatrix}.\end{aligned} \quad\quad\quad(29)

The Jacobian matrix has been written out explicitly as scalars because we are now switching to an abusive notation using matrices of vector elements. Our Jacobian, a matrix of scalars happened to have a nice compact representation in column vector form, but we cannot use this when multiplying out with our matrix elements (or perhaps could if we invented more conventions, but lets avoid that). Having written it out in full we see that we recover our original compact Jacobian representation, and have just

\begin{aligned}\boldsymbol{\nabla} = \begin{bmatrix}\hat{\mathbf{r}} & \frac{1}{{r}} \frac{\partial {\hat{\mathbf{r}}}}{\partial {\theta}} & \frac{1}{{r \sin^2\theta}} \frac{\partial {\hat{\mathbf{r}}}}{\partial {\phi}} \end{bmatrix}\begin{bmatrix}\partial_r \\ \partial_\theta \\ \partial_\phi\end{bmatrix}.\end{aligned} \quad\quad\quad(30)

Expanding this last product we have the gradient in its spherical polar representation

\begin{aligned}\boldsymbol{\nabla} = \begin{bmatrix}\hat{\mathbf{r}} \frac{\partial {}}{\partial {r}} + \frac{1}{{r}} \frac{\partial {\hat{\mathbf{r}}}}{\partial {\theta}} \frac{\partial {}}{\partial {\theta}} + \frac{1}{{r \sin\theta}} \frac{1}{{\sin\theta}} \frac{\partial {\hat{\mathbf{r}}}}{\partial {\phi}} \frac{\partial {}}{\partial {\phi}}\end{bmatrix}.\end{aligned} \quad\quad\quad(31)

With the labels

\begin{aligned}\hat{\boldsymbol{\theta}} &= \frac{\partial {\hat{\mathbf{r}}}}{\partial {\theta}} \\ \hat{\boldsymbol{\phi}} &= \frac{1}{{\sin\theta}} \frac{\partial {\hat{\mathbf{r}}}}{\partial {\phi}},\end{aligned} \quad\quad\quad(32)

(having confirmed that these are unit vectors), we have the final result for the gradient in this representation

\begin{aligned}\boldsymbol{\nabla} = \hat{\mathbf{r}} \frac{\partial {}}{\partial {r}} + \frac{1}{{r}} \hat{\boldsymbol{\theta}} \frac{\partial {}}{\partial {\theta}} + \frac{1}{{r \sin\theta}} \hat{\boldsymbol{\phi}} \frac{\partial {}}{\partial {\phi}}.\end{aligned} \quad\quad\quad(34)

Here the matrix delimiters for the remaining one by one matrix term were also dropped.

# General expression for gradient in orthonormal frames.

Having done the computation for the spherical polar case, we get the result for any orthonormal frame for free. That is just

\begin{aligned}\boldsymbol{\nabla} = \sum_i (\boldsymbol{\nabla} q_i) \frac{\partial {}}{\partial {q_i}}.\end{aligned} \quad\quad\quad(35)

From each of the gradients we can factor out a unit vector in the direction of the gradient, and have an expression that structurally has the same form as 34. Writing $\hat{\mathbf{q}}_i = (\boldsymbol{\nabla} q_i)/{\left\lvert{\boldsymbol{\nabla} q_i}\right\rvert}$, this is

\begin{aligned}\boldsymbol{\nabla} = \sum_i {\left\lvert{\boldsymbol{\nabla} q_i}\right\rvert} \hat{\mathbf{q}}_i \frac{\partial {}}{\partial {q_i}}.\end{aligned} \quad\quad\quad(36)

These individual direction gradients are not necessarily easy to compute. The procedures outlined in [1] are a more effective way of dealing with this general computational task. However, if we want, we can at proceed this dumb obvious way and be able to get the desired result knowing only how to apply the chain rule, and the Cartesian definition of the gradient.

# References

[1] F.W. Byron and R.W. Fuller. Mathematics of Classical and Quantum Physics. Dover Publications, 1992.