Addendum to Symmetry of the Riemann-Christoffel Tensor

This is an addendum to the previous post on the symmetries of the Riemann-Christoffel tensor. I received a comment that I should be able to make the arguments purely intrinsicly and without relying on the surface basis which is defined extrinsicly. This is a great comment and I wish that I had thought of it myself. So, here we go.

I think the only part of my original argument that relied on extrinsic characteristics was on the proof of the symmetry with the first two and last two indices. That argument relied on the definition of the surface Christoffel symbols as the dot-product of the partial derivative of the basis with the basis as follows $$ \begin{equation} \eqalign{ \Gamma^\gamma_{\alpha\beta} &= \pard{\vec{S_\alpha}}{S^\beta} \cdot \vec{S^\gamma} \cr \Gamma_{\gamma,\alpha\beta} &= \pard{\vec{S_\alpha}}{S^\beta} \cdot \vec{S_\gamma} } \label{eq:extchristoffel} \end{equation} $$ To fix this we will use the intrinsic definition of the Christoffel symbol $$ \begin{equation} \eqalign{ \Gamma_{\gamma,\alpha\beta} &= \frac{1}{2} \left( \pard{S_{\gamma\beta}}{S^\alpha} + \pard{S_{\gamma\alpha}}{S^\beta} - \pard{S_{\beta\alpha}}{S^\gamma} \right) } \label{eq:intchristoffel} \end{equation} $$ Here again is the definition of the Riemann-Christoffel tensor. The last two terms (the products of the Christoffels) were shown to have the required symmetry in the previous post so I won't repeat that here. I will only be looking at the red-boxed terms $$ \begin{equation} \eqalign{ R_{\gamma\delta\alpha\beta} &= \bbox[5px,border:2px solid red]{\pard{\Gamma_{\gamma,\beta\delta}}{S^\alpha} - \pard{\Gamma_{\gamma,\alpha\delta}}{S^\beta}} + \Gamma_{\epsilon,\gamma\beta}\Gamma^\epsilon_{\alpha\delta} - \Gamma_{\epsilon,\gamma\alpha}\Gamma^\epsilon_{\beta\delta} } \label{eq:rsx} \end{equation} $$ Applying the definition from \ref{eq:intchristoffel} to the highlighted terms gives us $$ \begin{equation} \eqalign{ \pard{\Gamma_{\gamma,\beta\delta}}{S^\alpha} &= \frac{1}{2} \pard{}{S^\alpha} \left( \pard{S_{\gamma\delta}}{S^\beta} + \pard{S_{\gamma\beta}}{S^\delta} - \pard{S_{\delta\beta}}{S^\gamma} \right) \cr \pard{\Gamma_{\gamma,\alpha\delta}}{S^\beta} &= \frac{1}{2} \pard{}{S^\beta} \left( \pard{S_{\gamma\delta}}{S^\alpha} + \pard{S_{\gamma\alpha}}{S^\delta} - \pard{S_{\delta\alpha}}{S^\gamma} \right) } \label{eq:twoterms} \end{equation} $$ Expressing the difference (while getting rid of the \( \frac{1}{2} \) factor) and applying the second differential operator $$ \begin{equation} \eqalign{ 2 \left( \pard{\Gamma_{\gamma,\beta\delta}}{S^\alpha} - \pard{\Gamma_{\gamma,\alpha\delta}}{S^\beta} \right) &= + \bbox[5px,border:2px solid red] { \pardd{S_{\gamma\delta}}{S^\alpha}{S^\beta} } + \bbox[5px,border:2px solid green]{ \pardd{S_{\gamma\beta}}{S^\alpha}{S^\delta} } - \bbox[5px,border:2px solid blue] { \pardd{S_{\delta\beta}}{S^\alpha}{S^\gamma} } - \bbox[5px,border:2px solid red] { \pardd{S_{\gamma\delta}}{S^\beta}{S^\alpha} } - \bbox[5px,border:2px solid blue] { \pardd{S_{\gamma\alpha}}{S^\beta}{S^\delta} } + \bbox[5px,border:2px solid green]{ \pardd{S_{\delta\alpha}}{S^\beta}{S^\gamma}} } \label{eq:expanded} \end{equation} $$ The two terms in the red boxes are the same because differentiation commutes (assuming continuity, differentiability, etc.). Since we subtract one from the other, they vanish into the mist.

The two terms in the blue boxes are self symmetric because if you make the required index swap you get the same thing you started with because of the symmetry of the metric tensor and commutative property of differentation.

The two terms in the green boxes are mutally symmetric. If you make the required index swap on one of them, you get the other one.

So along with the other items in the other post, we now have an intrinsic proof of the interchange symmetry of the first two and last two indices of the Riemann-Christoffel tensor. Here again is a link to the prior post

Thanks to 'Lemma' who I assume to be Prof. Grinfeld for the comment. I will indulge and give myself a Grinfeld number of #1.

Symmetry of the Riemann-Christoffel Tensor

In this episode, we delve into some of the symmetric and anti-symmetric properties of the Riemann-Christoffel tensor. This work corresponds to problem 245 in Chapter 12 of Introduction to Tensor Analysis and the Calculus of Moving Surfaces. The answer key says something like "we did this on the final exam." So my goal is to fill that hole.

Property #1

The first property is the easiest to show and follows directly from the definition and that is that if we switch \(\alpha\) and \(\beta\) (the 3rd and 4th indices) you get the negative value. This is anti-symmetry and it looks like this $$ R^\gamma_{\cdot\delta\beta\alpha} = -R^\gamma_{\cdot\delta\alpha\beta} $$ If you look at the definition, it is evident. $$ \begin{equation} (\nabla_\alpha\nabla_\beta - \nabla_\beta\nabla_\alpha)T^\gamma = R^\gamma_{\cdot\delta\alpha\beta}T^\delta \label{eq:commutator} \end{equation} $$ Switching the indices \(\alpha\) and \(\beta\) $$ \begin{equation} \eqalign { (\nabla_\beta\nabla_\alpha - \nabla_\alpha\nabla_\beta)T^\gamma &= R^\gamma_{\cdot\delta\beta\alpha}T^\delta \cr -(\nabla_\alpha\nabla_\beta - \nabla_\beta\nabla_\alpha)T^\gamma &= R^\gamma_{\cdot\delta\beta\alpha}T^\delta \cr -(\nabla_\alpha\nabla_\beta - \nabla_\beta\nabla_\alpha)T^\gamma &= -R^\gamma_{\cdot\delta\alpha\beta}T^\delta } \label{eq:commutator2} \end{equation} $$ Since this is true for an arbitrary tensor, \(T^\gamma\), we have $$ \begin{equation} R^\gamma_{\cdot\delta\beta\alpha} = -R^\gamma_{\cdot\delta\alpha\beta} \label{eq:sym1} \end{equation} $$

Property #2

The next symmetric property will require a little more finesse. Now, I don't claim that this is the easiest or most intuitive way to prove this but it does work and it is the one that I came up with. Besides, if you have a better way then it will still be good to see an alternative as that seems to broaden one's understanding. And perhaps you could be so kind as to give me a few clues in the comments.

Here we will show that the Riemann-Christoffel tensor with all indices lowered is symmetric if you swap the first two indices with the second two in order. In other words, \(\gamma\leftrightarrow\alpha\) and \(\delta\leftrightarrow\beta\). We start by going to the definition. $$ \begin{equation} \eqalign { R_{\gamma\delta\alpha\beta} &= \pard{\Gamma_{\gamma,\beta\delta}}{S^\alpha} - \pard{\Gamma_{\gamma,\alpha\delta}}{S^\beta} + \Gamma_{\epsilon,\gamma\beta}\Gamma^\epsilon_{\alpha\delta} - \Gamma_{\epsilon,\gamma\alpha}\Gamma^\epsilon_{\beta\delta} } \label{eq:rlower} \end{equation} $$ It's clear that the last term \(\Gamma_{\epsilon,\gamma\alpha}\Gamma^\epsilon_{\beta\delta}\) is symmetric with \(\{\alpha,\beta\}\leftrightarrow\{\gamma,\delta\}\) because the Christoffel symbols are symmetric in the last two indices. $$ \Gamma_{\epsilon,\gamma\alpha}\Gamma^\epsilon_{\beta\delta} = \Gamma_{\epsilon,\alpha\gamma}\Gamma^\epsilon_{\delta\beta} $$ One down, three to go. The third term requires a few additional steps. $$ \eqalign{ \Gamma_{\epsilon,\gamma\beta}\Gamma^\epsilon_{\alpha\delta} &= \left(S_{\omega\epsilon}\Gamma^\omega_{\gamma\beta}\right)\Gamma^\epsilon_{\alpha\delta} && \text{definition of lowered index} \\ &= \Gamma^\omega_{\gamma\beta} \left( S_{\omega\epsilon} \Gamma^\epsilon_{\alpha\delta} \right) && \text{commutativity/associativity} \\ &= \Gamma^\epsilon_{\gamma\beta} \Gamma_{\epsilon,\alpha\delta} && \text{lower the index of second term and rename the contracted index } \omega \rightarrow \epsilon \\ } $$ We see that the result is the same as the initial expression with the required indices swapped. Two down.

The last part of showing this symmetry is a bit more harrowing as we will need to consider the remaining two terms in tandem. We start by rewriting the first two terms of the definition of \(R\) in terms of the definitions of the Christoffel symbols of the surface. The definitions are $$ \begin{equation} \eqalign{ \Gamma^\gamma_{\alpha\beta} &= \pard{\vec{S_\alpha}}{S^\beta} \cdot \vec{S^\gamma} \cr \Gamma_{\gamma,\alpha\beta} &= \pard{\vec{S_\alpha}}{S^\beta} \cdot \vec{S_\gamma} } \label{eq:christoffel} \end{equation} $$

[EDIT: I've amended this proof with a fully intrinsic version of this portion of the proof. If this extrinsic argument gives you the willies, you can go here

Where \(\vec{S_\gamma}\) is the covariant basis of the surface. Rewriting the first two terms with these definitions and expanding via the product rule $$ \begin{equation} \eqalign { \pard{\Gamma_{\gamma,\beta\delta}}{S^\alpha} - \pard{\Gamma_{\gamma,\alpha\delta}}{S^\beta} &= \pard{}{S^\alpha} \left( \pard{\vec{S_\beta}}{S^\delta} \cdot \vec{S_\gamma} \right) - \pard{}{S^\beta} \left( \pard{\vec{S_\alpha}}{S^\delta} \cdot \vec{S_\gamma} \right) \cr &= \frac{\partial^2 \vec{S_\beta}}{\partial S^\alpha \partial S^\delta} \cdot \vec{S_\gamma} + \bbox[5px,border:2px solid red]{\pard{\vec{S_\beta}}{S^\delta} \cdot \pard{\vec{S_\gamma}}{S^\alpha}} - \frac{\partial^2 \vec{S_\alpha}}{\partial S^\beta \partial S^\delta} \cdot \vec{S_\gamma} - \bbox[5px,border:2px solid blue]{\pard{\vec{S_\alpha}}{S^\delta} \cdot \pard{\vec{S_\gamma}}{S^\beta}} } \label{eq:first2terms} \end{equation} $$ Exchanging indices \(\gamma\leftrightarrow\alpha\) and \(\delta\leftrightarrow\beta\) in the term boxed in blue yields the same two terms in reverse order. So, enough of that one.

The term boxed in red requires the realization that because of their definition the partial derivatives of the basis vectors have the property that \( \pard{\vec{S_\omega}}{S^\epsilon} = \pard{\vec{S_\epsilon}}{S^\omega} \). So, $$ \eqalign { \pard{\vec{S_\beta}}{S^\delta} \cdot \pard{\vec{S_\gamma}}{S^\alpha} &= \pard{\vec{S_\delta}}{S^\beta} \cdot \pard{\vec{S_\alpha}}{S^\gamma} } $$ Which is the same as swapping the required indices. Lastly, we examine the remaining two terms. "Factoring" out the derivative with respect to \( S^\delta \) gives us $$ \frac{\partial^2 \vec{S_\beta}}{\partial S^\alpha \partial S^\delta} \cdot \vec{S_\gamma} - \frac{\partial^2 \vec{S_\alpha}}{\partial S^\beta \partial S^\delta} \cdot \vec{S_\gamma} = \pard{}{S^\delta} \left( \pard{\vec{S^\beta}}{S^\alpha} - \pard{\vec{S^\alpha}}{S^\beta} \right) $$ By the same argument as the previous step, the difference inside the parentheses is identically zero. Now that all the terms are accounted for, we can state that the Riemann-Christoffel tensor is symmetric in \(\{\alpha,\beta\}\leftrightarrow\{\gamma,\delta\}\).

Property #3

Now that we have the above symmetries, showing the anti-symmetry of \( \gamma \leftrightarrow \delta \) is just 3 swaps away. $$ \eqalign { R_{\delta\gamma\beta\alpha} &= R_{\beta\alpha\delta\gamma} && \text{Property #2} \cr &= -R_{\beta\alpha\gamma\delta} && \text{Property #1} \cr &= -R_{\gamma\delta\alpha\beta} && \text{Property #2} \cr } $$

Values in the 2 Dimensional Riemann-Christoffel Tensor

The symmetries greatly restrict the degrees of freedom of the values in the tensor. Wikepedia tells me that the degrees of freedom from a "simple calculation" can be found to be $$ N = \frac{n^2(n^2 - 1)}{12} $$ In our case, \( n = 2 \) so we would expect one independent component. Here are all 16 values where \( x \) is the only value we can choose. $$ \eqalign { R_{1212} &= x \cr R_{1221} &= -x \cr R_{2112} &= -x \cr R_{2121} &= x \cr &\text{All other elements are 0} } $$

Lowering the Index on the Riemann-Christoffel Tensor

For the past few months I've been learning about Tensor Calculus. My original goal was to learn about Relativity but there was too much I didn't understand because of the use of tensors. I've always been intrigued by the word "tensors" but not having the foggiest idea what they were I had no ability to go forward.

So, I found a lecture series on YouTube called Tensor Calculus and the Calculus of Moving Surfaces taught by Professor Pavel Grinfeld. I was so intrigued that I bought the text book Introduction to Tensor Analysis and the Calculus of Moving Surfaces and proceeded to do the problems. I've just gotten to the chapter on curvature and figured I'd attempt solutions to some of the problems on-line. (I find the answer key somewhat lacking at times. Dr. Grinfeld has links on the YouTube channel for various resources.) My only quibble about the book is that there are a number of typographical errors. I hope it goes through a second printing.

The other thing I should note before I get started is that there is a lot of stuff on the web about differential geometry and tensors and the notation is all over the map. One of the pleasing things about Grinfeld's treatment is that the notation is fairly simple and consistent. So, I'll be using that notation because of this and also because it's the only one I know!

So, our problem today is to lower the index on the Riemann-Christoffel tensor. I don't pretend that this is the only or most affective or efficient way of doing this. It's the machinery that I employed because of the level of my understanding.

The Riemann-Christoffel tensor, \(R^\gamma_{\cdot\delta\alpha\beta}\) is the interesting components of the following expression; the commutator applied to a contravariant tensor. $$ \begin{equation} (\nabla_\alpha\nabla_\beta - \nabla_\beta\nabla_\alpha)T^\gamma = R^\gamma_{\cdot\delta\alpha\beta}T^\delta \label{eq:commutator} \end{equation} $$ Working out all the covariant derivatives you get the following famous expression for the Riemann-Christoffel tensor: $$ \begin{equation} R^\gamma_{\cdot\delta\alpha\beta} = \pard{\Gamma^\gamma_{\beta\delta}}{S^\alpha} - \pard{\Gamma^\gamma_{\alpha\delta}}{S^\beta} + \Gamma^\gamma_{\alpha\omega}\Gamma^\omega_{\beta\delta} - \Gamma^\gamma_{\beta\omega}\Gamma^\omega_{\alpha\delta} \label{eq:rup} \end{equation} $$ Lowering of the contravariant index, \(\gamma\), is achieved by multiplying by the metric tensor. $$ \begin{equation} \eqalign{ R_{\gamma\delta\alpha\beta} &= S_{\gamma\epsilon}R^\epsilon_{\cdot\delta\alpha\beta} \cr &= S_{\gamma\epsilon}\left(\pard{\Gamma^\epsilon_{\beta\delta}}{S^\alpha} - \pard{\Gamma^\epsilon_{\alpha\delta}}{S^\beta} + \Gamma^\epsilon_{\alpha\omega}\Gamma^\omega_{\beta\delta} - \Gamma^\epsilon_{\beta\omega}\Gamma^\omega_{\alpha\delta}\right) } \label{eq:rsx} \end{equation} $$ So, the interesting part of this is that we get accustomed to just dropping the indices because that is part of the beauty of tensor notation. The problem comes when you can't do that. The issue here is that multiplication of the metric tensor does not penetrate the partial derivative unscathed as it would for a covariant derivative. We have to be somewhat careful.

Here is what I did. Consider this expression which we can evaluate by the produce rule $$ \begin{equation} \pard{\Gamma_{\gamma,\beta\delta}}{S^\alpha} = \pard{\left(S_{\gamma\epsilon}\Gamma^\epsilon_{\beta\delta}\right)}{S^\alpha} = S_{\gamma\epsilon}\pard{\Gamma^\epsilon_{\beta\delta}}{S^\alpha} + \Gamma^\epsilon_{\beta\delta}\pard{S_{\gamma\epsilon}}{S^\alpha} \label{eq:lc2} \end{equation} $$ Note also that the partial derivative of the metric tensor can be found to be $$ \begin{equation} \pard{S_{\gamma\epsilon}}{S^\alpha} = \Gamma_{\gamma,\epsilon\alpha} + \Gamma_{\epsilon,\gamma\alpha} \label{eq:lc3} \end{equation} $$ Giving $$ \begin{equation} \eqalign{ S_{\gamma\epsilon}\pard{\Gamma^\epsilon_{\beta\delta}}{S^\alpha} &= \pard{\Gamma_{\gamma,\alpha\beta}}{S^\alpha} - \Gamma^\epsilon_{\beta\delta}\pard{S_{\gamma\epsilon}}{S^\alpha} \cr &= \pard{\Gamma_{\gamma,\beta\delta}}{S^\alpha} - \left(\Gamma_{\gamma,\epsilon\alpha} + \Gamma_{\epsilon,\gamma\alpha}\right)\Gamma^\epsilon_{\beta\delta} \cr &= \pard{\Gamma_{\gamma,\beta\delta}}{S^\alpha} - \Gamma_{\gamma,\epsilon\alpha}\Gamma^\epsilon_{\beta\delta} - \Gamma_{\epsilon,\gamma\alpha}\Gamma^\epsilon_{\beta\delta} } \label{eq:lc4} \end{equation} $$ If you look at equation \ref{eq:lc4}, it gives us the value of the first term of equation \ref{eq:rsx} after you multiply out the metric tensor. Doing the same thing to the second term gives us $$ \begin{equation} \eqalign{ S_{\gamma\epsilon}\pard{\Gamma^\epsilon_{\alpha\delta}}{S^\beta} &= \pard{\Gamma_{\gamma,\alpha\delta}}{S^\beta} - \Gamma_{\gamma,\epsilon\beta}\Gamma^\epsilon_{\alpha\delta} - \Gamma_{\epsilon,\gamma\beta}\Gamma^\epsilon_{\alpha\delta} } \label{eq:lc5} \end{equation} $$ The 3rd and 4th terms of equation \ref{eq:rsx} are just the Cristoffels with the index lowered. Putting everything together gives us $$ \begin{equation} \eqalign{ R_{\gamma\delta\alpha\beta} &= \pard{\Gamma_{\gamma,\beta\delta}}{S^\alpha} - \bbox[5px,border:2px solid red]{\Gamma_{\gamma,\epsilon\alpha}\Gamma^\epsilon_{\beta\delta}} - \Gamma_{\epsilon,\gamma\alpha}\Gamma^\epsilon_{\beta\delta} - \pard{\Gamma_{\gamma,\alpha\delta}}{S^\beta} + \bbox[5px,border:2px solid blue]{\Gamma_{\gamma,\epsilon\beta}\Gamma^\epsilon_{\alpha\delta}} + \Gamma_{\epsilon,\gamma\beta}\Gamma^\epsilon_{\alpha\delta} + \bbox[5px,border:2px solid red]{\Gamma_{\gamma,\alpha\omega}\Gamma^\omega_{\beta\delta}} - \bbox[5px,border:2px solid blue]{\Gamma_{\gamma,\beta\omega}\Gamma^\omega_{\alpha\delta}} } \label{eq:rl1} \end{equation} $$ Careful inspection will reveal that the red boxed terms are the same because the Cristoffel symbols are symmetric in the last two indices and the name of the contracted index doesn't matter. And the same goes for the terms in the blue boxes. It's interesting to note that the original Cristoffel products are disposed of and replaced with a different set. I don't know if there is some deeper meaning there.

The final result is (rearranging to get a more aesthetically pleasing result) $$ \begin{equation} \eqalign{ R_{\gamma\delta\alpha\beta} &= \pard{\Gamma_{\gamma,\beta\delta}}{S^\alpha} - \pard{\Gamma_{\gamma,\alpha\delta}}{S^\beta} + \Gamma_{\epsilon,\gamma\beta}\Gamma^\epsilon_{\alpha\delta} - \Gamma_{\epsilon,\gamma\alpha}\Gamma^\epsilon_{\beta\delta} } \label{eq:final} \end{equation} $$ More coming on the this topic.

Portland to New York Using only Gravity (Part 3)

Here we are at part 3 of a series of posts relating to an "underground brachistochone" where we wish to build a tunnel for a train propelled only by gravity from Portland, Oregon to New York, New York.

We left off in part 2 with the computation of the time to traverse the tunnel equal to $$ \Delta t = \pi \sqrt{\frac{R^2 - r_0^2}{gR}} $$ Of course, this requires us to know the depth of the tunnel but we were able to compute the time as a function of a fraction of the earth's radius, \(R\). What we want is a way to compute the depth of the tunnel, \(r\) as a function of the angle we're at or vice versa.

To start this process, we recall, from part 2 $$ \begin{equation} \theta' = \frac{1}{r}\sqrt{\frac{\frac{k^2g}{R}(R^2-r^2)}{r^2-\frac{k^2g}{R}(R^2-r^2)}} \label{eq:tp} \end{equation} $$ and recalling our definition of \(r_0\) stemming from solving for the conditions at the curve minimum $$ \begin{equation} \frac{k^2g}{R} = \frac{r_0^2}{R^2-r_0^2} \label{eq:r0} \end{equation} $$ Substuting \ref{eq:r0} into \ref{eq:tp} gives us, after a little bit of algebra, $$ \begin{equation} d\theta = \frac{1}{r}\sqrt{\frac {r_0^2(R^2-r^2)} {R^2(r^2-r_0^2)} }dr \label{eq:dt} \end{equation} $$ Moving some constants out of the way helps, sometimes, so lets do $$ \begin{equation} \frac{R}{r_0}d\theta = \frac{1}{r}\sqrt{\frac {R^2-r^2} {r^2-r_0^2} }dr \label{eq:dt2} \end{equation} $$ If you were to enter \ref{eq:dt2} into Wolfram-Alpha you would get a hot mess in return. Wolfram-Alpha is really, really smart but sometimes a little trickery is required to bend the integrand into something that it finds more recognizable. To that end, let $$ u^2 = \frac{r^2-r_0^2}{R^2-r^2} $$ This is kind of strange but making \(u\) equal to the upside-down-integrand kind of makes sense as the derivative of a square root is 1-over-the-square-root (multiplied by the derivative of the ... radicand ... it's been a while, is that right)? Solving for \(r^2\), $$ r^2 = \frac{u^2R^2 + r_0^2}{1+u^2} $$ Now we find \(dr\) but instead of taking the square root and differentiating, we will differentiate implicitly. $$ \eqalign { 2rdr &= {{2R^2u(1+u^2) - 2u(u^2R^2+r_0^2)} \over {(1+u^2)^2}}du \cr &= {{2u(R^2+r_0^2)} \over {(1+u^2)^2}}du \cr dr &= \frac{u}{r}{{(R^2+r_0^2)} \over {(1+u^2)^2}}du } $$ This will get a little messy but hang in there. First off, notice that $$ u=\sqrt{\frac{r^2 - r_0^2}{R^2-r^2}} $$ so $$ \frac{1}{u}=\sqrt{\frac{R^2-r^2}{r^2 - r_0^2}} $$ which is the most of the original integrand. Working our way through this $$ \eqalign { \frac{R}{r_0}\theta = \int\frac{1}{r}\sqrt{\frac {R^2-r^2} {r^2-r_0^2} }dr &= \int\frac{1}{r}\frac{1}{u}\frac{u}{r}\frac{R^2-r_0^2}{(1+u^2)^2}du \cr &= \int\frac{1}{r^2}\frac{R^2-r_0^2}{(1+u^2)^2}du \cr \text{(substituting for } r^2 \text{)} &= \int\frac{1+u^2}{(u^2R^2+r_0^2)}\frac{R^2-r_0^2}{(1+u^2)^2}du \cr &= \int\frac{R^2-r_0^2}{(u^2R^2+r_0^2)(1+u^2)}du \cr } $$ That last integrand looks like a job for partial fractions. $$ \eqalign { \frac{R^2-r_0^2}{(u^2R^2+r_0^2)(1+u^2)} &= \frac{A}{u^2R^2+r_0^2} + \frac{B}{1+u^2} \cr R^2-r_0^2 &= A(1+u^2) + B(u^2R^2+r_0^2) \cr &= u^2(A+BR^2) + A + Br_0^2 \cr \therefore A=R^2, B=-1 } $$ Putting this back into the integral $$ \begin{equation} \eqalign{ \frac{R}{r_0}\theta &= \int\frac{R^2-r_0^2}{(u^2R^2+r_0^2)(1+u^2)}du = R^2\int\frac{du}{u^2R^2+r_0^2} - \int\frac{du}{1+u^2} \cr &= \frac{R}{r_0}\tan^{-1}\left(\frac{Ru}{r_0}\right) - \tan^{-1}u + C \cr \theta &= \tan^{-1}\left(\frac{R}{r_0} \sqrt{\frac{r^2 - r_0^2}{R^2-r^2}} \right) - \frac{r_0}{R}\tan^{-1}\sqrt{\frac{r^2 - r_0^2}{R^2-r^2}} + C\cr } \end{equation} $$ Here we let \(\theta(r_0) = 0 \rightarrow C = 0\). At the surface, when \(r = R\), we can compute that $$ \theta = \frac{\pi}{2} - \frac{r_0}{R}\frac{\pi}{2} $$ $$ \begin{equation} \theta = \frac{\pi}{2} (1 - \frac{r_0}{R}) \label{eq:thetaR} \end{equation} $$ The angle \(\theta(R)\) in \ref{eq:thetaR} is the angle from where the curve is the deepest to the surface. So, given that the tunnel will travel from the surface down to the minimum then back up again, the total span, \(\Psi\), of the tunnel is $$ \Psi = 2\theta = \pi (1 - \frac{r_0}{R}) $$ And solving for \(r_0\) $$ \begin{equation} r_0 = R(1 - \frac{\Psi}{\pi}) \label{eq:r0x} \end{equation} $$ The earth's radius in kilometers is \(R=6371\text{km}\). The flight distance from PDX to JFK is \(d=3949\text{km}\). Since \(d=R\Psi\) the required angle spanned by this distance is \(\Psi = \frac{d}{R} = \frac{3949}{6371} = 0.62\). $$ \eqalign { \Psi &= 0.62 \cr r_0 &= R(1 - \frac{\Psi}{\pi}) = 6371(1 - \frac{0.62}{\pi}) = 5114\text{ km} \cr } $$ And from part 2 of this series, the time \(\Delta t\) is given by (remember to adjust the units!) $$ \eqalign { \Delta t &= \pi \sqrt{\frac{R^2 - r_0^2}{gR}} \cr &= \pi \sqrt{\frac{6371000^2 - 5114000^2}{(9.8)(6371000)}} \cr &\approx 25 \text{ minutes} } $$ We can also ask about the maximum speed attained. Since we've declared that the sum of the kinetic and potential energy are zero and that this is true for the entire traversal of the path, we can deduce that the kinetic energy is at a maximum where the potential energy is at a minimum. This happens at \(r=r_0\). We saw from part one of this series that $$ \eqalign { |v| &= \sqrt{\frac{g}{R}(R^2-r^2)} \cr &= \sqrt{\frac{9.8}{6371000}(6371000^2 - 5114000^2)} \cr &\approx 4713 \frac{\text{m}}{\text{s}} \cr &\approx 16970 \frac{\text{km}}{\text{hr}} } $$ So there you have it... How long is the tunnel? How many Gs are being generated at \(r=r_0\). I'll leave those for the reader.

Underground from Portland to New York (Part 2)

This post is part two of what will be a three part series on the derivation of the equations of an underground tunnel that minimizes travel time when only gravity is used for propulsion. We left our subterranean tunnel problem with the following equation we wished to minimize. $$ \eqalign{ I &= \int_a^b\frac{\sqrt{(dr)^2+(rd\theta)^2}}{\sqrt{\frac{g}{R}(R^2-r^2)}} } $$ Our first issue is figuring out what to do with \(\sqrt{(dr)^2+(rd\theta)^2}\) so that we can integrate it. We could choose to factor out either the \(dr\) or the \(d\theta\) and could solve the problem either way. We'll choose to factor out the \(dr\) because that leaves us an integrand that is not a function of the path variable (in this case \(\theta\)) which means that \(\pard{F}{\theta} = 0\) and \(\pard{F}{\theta'} = constant\). So, we have the Euler-Lagrange equation for this problem (remember that \(r\) is our variable of integration): $$ \eqalign{ \frac{d}{dr}\pard{F}{\theta'} - \pard{F}{\theta} = 0 \cr \text{but}\enspace \pard{F}{\theta} = 0 \cr \text{therefore}\enspace \frac{d}{dr}\pard{F}{\theta'} = 0 \cr \text{so}\enspace \pard{F}{\theta'} = k } $$ $$ \pard{}{\theta'} \left(\frac{\sqrt{1+(r\theta')^2}}{\sqrt{\frac{g}{R}(R^2-r^2)}}\right) = k $$ Performing the differentiation we have $$ \frac{1}{\sqrt{\frac{g}{R}(R^2-r^2)}}\cdot\frac{r^2\theta'}{\sqrt{1+(r\theta')^2}} = k $$ Squaring both sides and solving for \((r\theta')^2\) first $$ \frac{1}{\frac{g}{R}(R^2-r^2)}\cdot\frac{(r^2\theta')^2}{1+(r\theta')^2} = k^2 $$ $$ \begin{equation} \frac{\frac{k^2g}{R}(R^2-r^2)}{r^2-\frac{k^2g}{R}(R^2-r^2)} = (r\theta')^2 \label{eq:rt2} \end{equation} $$ We solve for \((r\theta')^2\) and mark that equation because in the end, our original integral for the minimum time is a function of \((r\theta')^2\) and we don't want to have to recompute it. Going the final step to solve for \(\theta'\) we have $$ \begin{equation} \theta' = \frac{1}{r}\sqrt{\frac{\frac{k^2g}{R}(R^2-r^2)}{r^2-\frac{k^2g}{R}(R^2-r^2)}} \label{eq:tp} \end{equation} $$ At this point we pause to consider the variable and limits of integration. If you look at the figure below, we can see that to traverse the curve, \(r\) goes from \(R\) to \(R\). That doesn't help us much because the definite integral would end up being 0.

We do notice, however, that the curve must be symmetric (the path shouldn't matter on which side you start) and so the deepest part of the curve is in the middle at \(r_0\). At this value of \(r=r_0\), we can see that \(\frac{dr}{d\theta} = 0\) since \(r(\theta)\) is mimimized. Conversely, this means that \(\theta' = \frac{d\theta}{dr} \rightarrow \infty\). This happens when the denominator of equation \ref{eq:tp} becomes zero. So, solving for the constants in this case when \(r=r_0\) $$ r_0^2 - \frac{k^2g}{R}(R^2-r_0^2) = 0 $$ yields $$ \frac{k^2g}{R} = \frac{r_0^2}{R^2-r_0^2} $$ Before we find the path, we'll find the time it takes to traverse the curve. Of course, this time will be a function of \(r_0\) so we will still have some work to do. Rewriting equation \ref{eq:rt2} with the above substitution and then cleaning up some of the fractions yields $$ (r\theta')^2 = \frac{\frac{r_0^2}{R^2-r_0^2}(R^2-r^2)}{r^2-\frac{r_0^2}{R^2-r_0^2}(R^2-r^2)} $$ yields $$ \begin{equation} (r\theta')^2 = \frac{r_0^2(R^2-r^2)}{r^2(R^2-r_0^2)-r_0^2(R^2-r^2)} \label{eq:rt2r0} \end{equation} $$ Now, going back to our original functional we factor out the \(dr\) and put in the actual limits of integration, multiplying the whole thing by 2 because we're only going halfway in the integral. $$ \Delta t = 2\int_{R}^{r_0}\frac{\sqrt{1+(r\theta')^2}}{\sqrt{\frac{g}{R}(R^2-r^2)}}dr $$ Substituting \ref{eq:rt2r0} into the above and moving some constants over to the left hand side $$ \frac{1}{2}\sqrt{\frac{g}{R}}\Delta t = \int_{R}^{r_0}\frac{1}{\sqrt{(R^2-r^2)}}\cdot\frac{r\sqrt{R^2-r_0^2}}{\sqrt{r^2(R^2-r_0^2) - r_0^2(R^2-r^2)}}dr $$ And after a little more algebra, factoring out the \((R^2-r_0^2)\) $$ \frac{1}{2}\sqrt{\frac{g}{R}}\Delta t = \int_{R}^{r_0}\frac{r}{\sqrt{(R^2-r^2)(r^2 - \frac{r_0^2}{R^2-r_0^2}(R^2-r^2))}}dr $$ Sometimes at points like this it helps to look at the units to make sure we at least have a sanity check. On the left hand side we have \(\sqrt{\frac{m}{s^2}\cdot\frac{1}{m}}\cdot s\) which is unitless. On the right, we have \(\frac{m\cdot m}{\sqrt{m^4}}\) which is also unitless. So we at least have that going for us. We're almost home for this part. We now have a relatively straightforward integration. We let \(u=R^2-r^2\), \(du = -2rdr\) and \(r^2 = R^2 - u\). These substitutions result in $$ \frac{1}{2}\sqrt{\frac{g}{R}}\Delta t = -\frac{1}{2}\int_{R}^{r_0}\frac{1}{\sqrt{u((R^2-u) - \frac{r_0^2}{R^2-r_0^2}u)}}dr $$ Simplifying the denominator $$ \frac{1}{2}\sqrt{\frac{g}{R}}\Delta t = -\frac{1}{2}\int_{R}^{r_0}\frac{1}{\sqrt{u(R^2-(1+\frac{r_0^2}{R^2-r_0^2})u)}}dr $$ A little more $$ \frac{1}{2}\sqrt{\frac{g}{R}}\Delta t = -\frac{1}{2}\int_{R}^{r_0}\frac{1}{\sqrt{u(R^2-\frac{R^2}{R^2-r_0^2}u)}}dr $$ A little more algebra by factoring out an \(R^2\) from the to get the integral into a more recognizable form $$ \frac{1}{2}\sqrt{\frac{gR^2}{R}}\Delta t = -\frac{1}{2}\int_{R}^{r_0}\frac{1}{\sqrt{u(1 - \frac{1}{R^2-r_0^2}u)}}dr $$ Simplifying the left a little and with a little help from Wolfram-Alpha we have (putting back the \(u=R^2-r^2\)) $$ \eqalign { \frac{1}{2}\sqrt{gR}\Delta t &= 2\sqrt{(R^2-r_0^2)}\sin^{-1}\left(\sqrt{\frac{R^2-r^2}{R^2-r_0^2}}\right) \Bigg \bracevert_R^{r_0} \cr \frac{1}{2}\sqrt{gR}\Delta t &= 2\sqrt{(R^2-r_0^2)}(\frac{\pi}{2} - 0) \cr \frac{1}{2}\sqrt{gR}\Delta t &= \pi\sqrt{R^2-r_0^2} } $$ And finally solving for \(\Delta t\) $$ \Delta t = 2 \pi \sqrt{\frac{R^2 - r_0^2}{gR}} $$ Here, we can ask the question of how long will it take to go from one side of the earth to the other, straight through the core. At that point, \(r_0 = 0\) so \(\Delta t = 2\pi\sqrt{\frac{R}{g}} \). With \(g=9.8\frac{\text{m}}{\text{s}^2}\) and the radius of the earth \(R=6.4\cdot10^6 \text{m}\) $$ \Delta t = 2\pi\sqrt{\frac{6.4\cdot 10^6}{9.8}} \approx 5078 \text{s} \approx 85 \text{min} $$ Here is a graph of the transit time in minutes versus the depth as a fraction of \(R\).

Solving for the path is coming in part 3.

Underground from Portland to New York (Part 1)

The Brachistochrone

The brachistochrone is a classic problem put forth by Johann Bernoulli and solved by some of the mathematical heavyweights of the day. The basic premise of the problem is that given two points connected by a curve, find the shape of the curve such that a particle moving only under the influence of gravity would travel along the curve in the minimum time. Somewhat surprisingly, this curve is not a straight line. Even Galileo figured this out in 1638 though mistakenly thought the curve was a circular arc. It turns out that the curve is a cycloid but that's not what we're going to prove here.

Going Underground

Lets say we want to build a transportation tunnel from Portland, Oregon to New York, New York. We've invented a frictionless surface and we can evacuate the tunnel so there won't be any air resistance. We want to go completely green and only use gravity for travel. Oh, and to add to our perfectly realistic assumptions, we live on an earth of uniform density and of perfectly spherical shape.

What is the shape of the tunnel that will get us across in the shortest amount of time?

In similar fashion to solving the Brachistochrone problem, we must minimize the functional $$ I=\int_{t0}^{t1}dt = \int_a^b \frac{ds}{v} $$

To get started we'll employ some basic physics. Because there is no energy added to the system, the sum of kinetic and potential energies will be a constant throughout the trip. The kinetic energy, \(T\), is familiar and we can write it directly as $$ T = \frac{1}{2}m_Tv^2 $$ where \(m_T\) is the mass of the train. The potential energy, \(V\), is not quite as straightforward as it is not the familiar \(V=mgh\) we've seen in our entry level physics class. To derive it we first need to find the forces involved. The force of gravity is $$ \vec{F} = -G\frac{m_Em_T}{r^2} \uvec{r} $$ where \(m_E\) is the effective mass of the earth when we are somewhere inside of it at a distance \(r\) from the center. This mass is all the mass that is located less than or equal to a distance \(r\) from the center. So $$ m_E = \frac{4}{3}\rho\pi r^3 $$ where \(\rho\) is the (uniform!) density of the earth. Plugging that back into our force equation gives $$ \vec{F} = -G\frac{\left(\frac{4}{3}\rho\pi r^3\right)m_T}{r^2} \uvec{r} = -\frac{4}{3}G\rho\pi m_T r \uvec{r} $$ (Note the negative sign since the force is directed opposite the radial direction.) Since we know the acceleration due to gravity is \(g = 9.8 \frac{m}{s^2} \) at the surface, \(r = R\) we make the substutition of $$ \eqalign{ g &= \frac{4}{3}G\pi\rho R \cr \vec{F} &= -\frac{m_T g}{R}r\uvec{r} } $$ As an aside, here it is interesting to note that inside the earth, the force actually increases linearly with the distance \(r\). This is because the mass increases as \(r^3\) even though the gravity is an inverse-square law.

Now to compute the integral for potential energy, we will choose our reference (\(V=0\)) point to be \(r = R\). The line integral is along the radius to a distance \(r\) from the center. $$ V = -\int_R^r -\frac{m_T g}{R} x dx = \frac{m_Tg}{2R}(r^2 - R^2) $$ Since the total energy \(T+V=0\) at the start and no energy is added to the system we can set the total energy to 0 and solve for the speed as a function of \(r\). $$ \eqalign{ T+V &= \frac{1}{2}m_Tv^2 + \frac{m_Tg}{2R}(r^2 - R^2) = 0 \cr v &= \sqrt{\frac{g}{R}(R^2-r^2)} } $$ We are now at a point where we can write the functional we wish to minimize as $$ \eqalign{ I &= \int_{t0}^{t1}dt = \int_a^b \frac{ds}{v} \cr &= \int_a^b\frac{\sqrt{(dr)^2+(rd\theta)^2}}{\sqrt{\frac{g}{R}(R^2-r^2)}} } $$ To be continued in part deux.

Derivation of Divergence in Cylindrical Coordinates

$$ \newcommand{\pard}[2]{\frac{\partial #1}{\partial #2}} $$ $$ \newcommand{\vec}[1]{\boldsymbol{#1}} $$ $$ \newcommand{\unitv}[1]{\boldsymbol{\hat{#1}}} $$ The divergence of a vector field in Cartesian coordinates is fairly straight forward: $$ \nabla\cdot \vec{F} = \pard{F_x}{x} + \pard{F_y}{y} + \pard{F_z}{z} $$ But before we undertake this derivation, it will help to try and understand the meaning of divergence. It is tempting to just take the partial derivatives in whatever coordinate system you're using and add them together. But alas, this doesn't quite work.

First, let us consider flux.

Flux is the amount of something passing through a surface multiplied by the area of that surface. It can be a fluid, gas, electric field, or anything that can be described by a vector field. In the picture below, there are F liters/minute of water flowing in the circular pipe.

If we were to consider an imaginary surface slicing perpendicularly through the pipe (surface A), the flux through that circular surface would be F liters/minute multiplied by the area of the circle. It is important to note that that if we were to consider a surface slicing through at a bias (like surface B), the flux would be the same. Even though the elliptical area is larger than the circular one the amount of flux per unit area is less because the flux is no longer perpendicular to the surface. Flux through an area is then defined as $$ \Psi = \vec{F}\cdot\vec{A} = \|\vec{F}\|\|\vec{A}\|\cos\theta $$ where the vector area \( \vec{A} \) is defined to be normal to the surface. For a surface that is open, the choice of which of the two normal vectors we mean is arbitrary (we just have to be consistent). For a closed surface, the normal is defined to point to the outside of the surface. The angle \(\theta\) is the angle between the area-normal and the vector \(\vec{F}\). You can see as the angle between the flow and the surface gets closer to \(\frac{\pi}{2}\) that the flux gets smaller and smaller until it vanishes completely when the two are perpendicular.

Now, consider the volume defined by surfaces A, B and the walls of the pipe in the figure above. We can define the total flux through the volume as $$ \Psi_{total} = \Psi_{pipe} + \Psi_A + \Psi_{B} $$ Unless we have a leaky pipe, \(\Psi_{pipe} = 0\). Because we define the surface normal of a closed surface to point outward, the flux through surface A is negative (the flow and the normal are in opposite directions, so the dot-product is negative). The flow through B is equal in magnitude to that passing through surface B but positive. Hence, \(\Psi_{total} = 0 \). This makes intuitive sense when dealing with conservation of mass as otherwise there would be an accumulation of mass inside the pipe (\(\Psi_{total} < 0\)) or loss of mass (\(\Psi_{total} > 0\)).

Divergence is the amount of flux leaving a closed surface as the volume of the closed surface approaces zero. So, we now look at the general sitution of the pipe above in cylindrical coordinates. We consider a small enough closed volume that we can consider it a cube. In order to compute the total flux, we have to do 6 integrations as there are 6 sides to a cube. But we will group them in pairs, one pair for each direction \(\rho, \phi, z\).

The vector field in cylindrical coordinates is given by $$ \vec{F} = F_{\rho} \unitv{\rho} + F_{\phi} \unitv{\phi} + F_{z} \unitv{z} $$ We place our tiny cube-like volume in the field \(\vec{F}\) at some point \(\vec{P}=(\rho,\phi,z)\).

First we consider the flux passing through the surfaces normal to the \(z\) direction (the top and bottom faces of the cube). Since the distances involved are small, we can approximate the field passing through these as $$ F_{top} = F_z(\rho,\phi,z + \Delta z) \approx F_z(\rho,\phi,z) + \pard{F_z}{z}(\rho,\phi,z)\frac{\Delta z}{2} $$ It should be noted that the other components of the field don't play a part here because they are parallel to these surfaces so nothing going in the, say, \(\rho\) direction leaves via the \(z\) direction. Similarly for the bottom surface, $$ F_{bottom} = F_z(\rho,\phi,z - \Delta z) \approx F_z(\rho,\phi,z) - \pard{F_z}{z}(\rho,\phi,z)\frac{\Delta z}{2} $$ The areas of these surfaces are equal and independent of \(z\). $$ A_z = \rho\Delta\phi\Delta\rho $$ Therefore, the flux through these z-surfaces is (remembering that the flux through the bottom surface is into the volume and so is negative) $$ \Psi_z = F_{top}A_z - F_{bottom}A_z = (F_{top} - F_{bottom})A_z $$ $$ \Psi_z = (F_z + \pard{z}{F_z}\frac{\Delta z}{2} - (F_z - \pard{F_z}{z}\frac{\Delta z}{2}))\rho\Delta\phi\Delta\rho $$ $$ \Psi_z = \pard{F_z}{z}\rho\Delta \phi \Delta\rho \Delta z $$ Recognize that the term \( \rho\Delta \phi \Delta\rho \Delta z \) is the differential volume in cylindrical coordinates, so we can write the previous equation as $$ \begin{equation} \label{eq:fluxz} \Psi_z = \pard{F_z}{z}\Delta V \end{equation} $$ Where \(\Delta V\) is the differential volume element. I should note that if you were to do this derivation in Cartesian coordinates, the x and y component derivations would be identical (except for the subscripts). But in cylindrical coordinates, things get a little more interesting from here on out. Looking now at the \(\phi\) direction we consider the front and back faces of the volume. (I will omit the function arguments to the partial derivatives in the service of brevity with the understanding that we are talking about the point \(P=(\rho,\phi,z)\)). $$ F_{back} = F_\phi(\phi + \Delta \phi) \approx F_\phi + \pard{F_\phi}{\phi}\frac{\Delta \phi}{2} $$ $$ F_{front} = F_\phi(\phi - \Delta \phi) \approx F_\phi - \pard{F_\phi}{\phi}\frac{\Delta \phi}{2} $$ $$ A_\phi = \Delta\rho\Delta z $$ Because the positive \(\phi\) direction is toward the back of the volume, $$ \Psi_{\phi} = (F_{back} - F_{front})A_{\phi} $$ $$ \Psi_{\phi} = (F_\phi + \pard{F_\phi}{\phi}\frac{\Delta \phi}{2} - (F_\phi - \pard{F_\phi}{\phi}\frac{\Delta \phi}{2}))\Delta\rho\Delta z $$ $$ \Psi_{\phi} = \pard{F_\phi}{\phi}\Delta\phi\Delta\rho\Delta z $$ The term \(\Delta\phi\Delta\rho\Delta z\) is almost the differential volume. And maybe because we know where we are headed, we bend the equation to our will $$ \Psi_{\phi} = \pard{F_\phi}{\phi}\Delta\phi\Delta\rho\Delta z = \pard{F_\phi}{\phi}\frac{\rho}{\rho}\Delta\phi\Delta\rho\Delta z $$ $$ \begin{equation} \label{eq:fluxphi} \Psi_\phi = \frac{1}{\rho}\pard{F_\phi}{\phi}\Delta V \end{equation} $$ Now we try to handle the \(\rho\) component. We need to be careful here. But first things first. The vector field is defined similarly to previous components. $$ F_{right} = F_\rho(\rho + \Delta \rho) \approx F_\rho + \pard{F_\rho}{\rho}\frac{\Delta \rho}{2} $$ $$ F_{left} = F_\phi(\rho - \Delta \rho) \approx F_\rho - \pard{F_\rho}{\rho}\frac{\Delta \rho}{2} $$ Here is where it gets tricky. The area definitions are different because the two surfaces are not equal in size and so must be treated differently. $$ A_{right} = (\rho+\frac{\Delta\rho}{2})\Delta\phi\Delta z $$ $$ A_{left} = (\rho-\frac{\Delta\rho}{2})\Delta\phi\Delta z $$ First computing the flux on the right side. $$ \Psi_{right} = F_{right}A_{right}=(F_\rho + \pard{F_\rho}{\rho}\frac{\Delta \rho}{2})((\rho+\frac{\Delta\rho}{2})\Delta\phi\Delta z) $$ $$ \Psi_{right} = (F_\rho + \pard{F_\rho}{\rho}\frac{\Delta \rho}{2})(\rho\Delta\phi\Delta z + \frac{1}{2}\Delta\rho\Delta\phi\Delta z) $$ $$ \Psi_{right} = F_\rho\rho\Delta\phi\Delta z+\frac{1}{2}F_\rho\Delta\rho\Delta\phi\Delta z + \frac{1}{2}\pard{F_\rho}{\rho}\rho\Delta\rho\Delta\phi\Delta z + \frac{1}{4}\pard{F_\rho}{\rho}\Delta\rho^2\Delta\phi\Delta z $$ Similarly on the left hand side $$ \Psi_{left} = F_{left}A_{left}=(F_\rho - \pard{F_\rho}{\rho}\frac{\Delta \rho}{2})((\rho-\frac{\Delta\rho}{2})\Delta\phi\Delta z) $$ $$ \Psi_{left} = (F_\rho - \pard{F_\rho}{\rho}\frac{\Delta \rho}{2})(\rho\Delta\phi\Delta z - \frac{1}{2}\Delta\rho\Delta\phi\Delta z) $$ $$ \Psi_{left} = F_\rho\rho\Delta\phi\Delta z-\frac{1}{2}F_\rho\Delta\rho\Delta\phi\Delta z - \frac{1}{2}\pard{F_\rho}{\rho}\rho\Delta\rho\Delta\phi\Delta z + \frac{1}{4}\pard{F_\rho}{\rho}\Delta\rho^2\Delta\phi\Delta z $$ Computing \(\Psi_\rho\) (finally, after all that algebra) by subtracting and cancelling terms $$ \Psi_\rho = \Psi_{right} - \Psi_{left} $$ $$ \Psi_\rho = F_\rho\Delta\rho\Delta\phi\Delta z + \pard{F_\rho}{\rho}\rho\Delta\rho\Delta\phi\Delta z $$ We make the same modification to the first term that we did with \(\Psi_\phi\) and substituting \(\Delta V\) $$ \begin{equation} \label{eq:fluxrho} \Psi_\rho = \frac{1}{\rho}F_\rho\Delta V + \pard{F_\rho}{\rho}\Delta V \end{equation} $$ This last equation is correct but we normally don't see it written this way, so we factor out a \(\frac{1}{\rho}\) and invoke the reverse chain rule. $$ \Psi_\rho = \frac{1}{\rho}(F_\rho + \rho\pard{F_\rho}{\rho})\Delta V = \frac{1}{\rho}\pard{}{\rho}(\rho F_\rho)\Delta V $$ And so the total flux through our tiny volume is $$ \Psi_{total} = \Psi_\rho + \Psi_\phi + \Psi_z = \left(\frac{1}{\rho}\pard{}{\rho}(\rho F_\rho) + \frac{1}{\rho}\pard{F_\phi}{\phi} + \pard{F_z}{z}\right)\Delta V $$ The divergence is the ratio of the flux through a closed surface as the volume of that closed surface approaces 0. $$ \nabla\cdot\vec{F} = \lim\limits_{\Delta V \to 0} \frac{\Psi_{total}}{\Delta V} = \frac{1}{\rho}\pard{}{\rho}(\rho F_\rho) + \frac{1}{\rho}\pard{F_\phi}{\phi} + \pard{F_z}{z} $$

The derivation of divergence in spherical coordinates procedes along similar lines as the \(\phi\) and \(\rho\) terms. It's a little more complicated because of the added \(\sin\theta\) terms.