Peeter Joot's (OLD) Blog.

Math, physics, perl, and programming obscurity.

Posts Tagged ‘cylindrical coordinates’

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

Posted by peeterjoot on March 3, 2013

In A compilation of notes, so far, for ‘PHY452H1S Basic Statistical Mechanics’ I posted a link this compilation of statistical mechanics course notes.

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

Posted in Math and Physics Learning. | Tagged: , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , | 1 Comment »

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]


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{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)


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)


\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{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)


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.


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

Posted in Math and Physics Learning. | Tagged: , , , , , | Leave a Comment »

Couette flow of a viscous incompressible fluid (visualizing it).

Posted by peeterjoot on April 12, 2012

[Click here for a PDF of this post(and previous).]

In Couette flow of a viscous incompressible fluid, I looked at solving this problem. As a followup, here are some visualizations of the velocity field as the outer cylinder angular velocity is changed.

The Mathematica workbook that was used to generate the figures is also available for more direct play, and includes slider controls for the cylinder radii and angular velocities.

Posted in Math and Physics Learning. | Tagged: , , , , , , | Leave a Comment »

Couette flow of a viscous incompressible fluid.

Posted by peeterjoot on April 9, 2012

[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.)]


A problem from this years phy1530 problem set 2 that appears appropriate for phy454 exam prep.


Consider the steady flow between two long cylinders of radii R_1 and R_2, R_1 > R_1, rotating about their axes with angular velocities \Omega_1, \Omega_2. Look for a solution of the form, where \hat{\boldsymbol{\phi}} is a unit vector along the azimuthal direction:


\begin{aligned}\mathbf{u} = v(r) \hat{\boldsymbol{\phi}}\end{aligned} \hspace{\stretch{1}}(2.1a)

\begin{aligned}p = p(r).\end{aligned} \hspace{\stretch{1}}(2.1b)


  • Write out the Navier-Stokes equations and find differential equations for v(r) and p(r). You should find that these equations have relatively simple solutions, i.e.,

    \begin{aligned}v(r) = a r + \frac{b}{r}.\end{aligned} \hspace{\stretch{1}}(2.2)

  • Fix gthe constants a and b from the boundary conditions. Determine the pressure p(r).
  • Compute the friction forces that the fluid exerts on the cylinders, and compute the torque on each cylinder. Show that the total torque on the fluid is zero (as must be the case).

This is also a problem that I recall was outlined in section 2 from [1]. Some of the instabilities that are mentioned in the text are nicely illustrated in [2].

We illustrate our system in figure (1).

Figure 1: Coutette flow configuration


Solution: Part 1. Navier-Stokes and resulting differential equations.

Navier-Stokes for steady state incompressible flow has the form


\begin{aligned}(\mathbf{u} \cdot \boldsymbol{\nabla}) \mathbf{u} = -\frac{1}{{\rho}} \boldsymbol{\nabla} p + \nu \boldsymbol{\nabla}^2 \mathbf{u}\end{aligned} \hspace{\stretch{1}}(3.3a)

\begin{aligned}\boldsymbol{\nabla} \cdot \mathbf{u} = 0.\end{aligned} \hspace{\stretch{1}}(3.3b)


where the gradient has the form

\begin{aligned}\boldsymbol{\nabla} = \hat{\mathbf{r}} \partial_r + \frac{\hat{\boldsymbol{\phi}}}{r} \partial_\phi.\end{aligned} \hspace{\stretch{1}}(3.4)

Let’s first verify that the incompressible condition 3.3b is satisfied for the presumed form of the solution we seek. We have

\begin{aligned}\boldsymbol{\nabla} \cdot \mathbf{u} &=\left( \hat{\mathbf{r}} \partial_r + \frac{\hat{\boldsymbol{\phi}}}{r} \partial_\phi \right) \cdot (v(r) \hat{\boldsymbol{\phi}}(\phi) ) \\ &=(\hat{\mathbf{r}} \cdot \hat{\boldsymbol{\phi}}) v' + \frac{\hat{\boldsymbol{\phi}}^2}{r} \partial_\phi v(r)+ \frac{v(r) \hat{\boldsymbol{\phi}}}{r} \cdot \partial_\phi \hat{\boldsymbol{\phi}} \\ &= \frac{v(r) \hat{\boldsymbol{\phi}}}{r} \cdot (-\hat{\mathbf{r}}) \\ &= 0\end{aligned}

Good. Now let’s write out the terms of the momentum conservation equation 3.3a. We’ve got

\begin{aligned}(\mathbf{u} \cdot \boldsymbol{\nabla}) \mathbf{u}&=\frac{ v}{r} \partial_\phi ( v \hat{\boldsymbol{\phi}} ) \\ &=-\frac{ v^2 \hat{\mathbf{r}}}{r},\end{aligned}


\begin{aligned}-\frac{1}{{\rho}} \boldsymbol{\nabla} p&=-\frac{1}{{\rho}} \left( \hat{\mathbf{r}} \partial_r + \frac{\hat{\boldsymbol{\phi}}}{r} \partial_\phi \right) p(r) \\ &=-\frac{\hat{\mathbf{r}} p'}{\rho},\end{aligned}


\begin{aligned}\nu \boldsymbol{\nabla}^2 \mathbf{u}&=\nu \left( \hat{\mathbf{r}} \partial_r + \frac{\hat{\boldsymbol{\phi}}}{r} \partial_\phi \right) \cdot\left( \hat{\mathbf{r}} \partial_r + \frac{\hat{\boldsymbol{\phi}}}{r} \partial_\phi \right) (v(r) \hat{\boldsymbol{\phi}}(\phi)) \\ &=\nu \left( \partial_{rr} + \frac{1}{{r^2}} \partial_{\phi\phi}+ \frac{\hat{\boldsymbol{\phi}}}{r} \partial_\phi \cdot (\hat{\mathbf{r}} \partial_r)\right)(v(r) \hat{\boldsymbol{\phi}}(\phi)) \\ &=\nu \left( \partial_{rr} + \frac{1}{{r^2}} \partial_{\phi\phi}+ \frac{1}{r} \partial_r\right)(v(r) \hat{\boldsymbol{\phi}}(\phi)) \\ &=\nu \left( \frac{1}{{r}} \partial_{r} (r \partial_r) + \frac{1}{{r^2}} \partial_{\phi\phi}\right)(v(r) \hat{\boldsymbol{\phi}}(\phi)) \\ &=\nu \left( \frac{1}{{r}} (r v')' - \frac{v}{r^2} \right)\hat{\boldsymbol{\phi}}\end{aligned}

Equating \hat{\mathbf{r}} and \hat{\boldsymbol{\phi}} components we have two equations to solve


\begin{aligned}r (r v')' - v = 0\end{aligned} \hspace{\stretch{1}}(3.5a)

\begin{aligned}p' = \frac{\rho v^2}{r}.\end{aligned} \hspace{\stretch{1}}(3.5b)


Expanding out our velocity equation we have

\begin{aligned}r^2 v'' + r v' - v = 0,\end{aligned} \hspace{\stretch{1}}(3.6)

for which we’ve been told to expect that 2.2 is a solution (and it has the two integration constants we require for a solution to a homogeneous equation of this form). Let’s verify that we’ve computed the correct differential equation for the problem by trying this solution

\begin{aligned}r^2 v'' + r v' - v &=r^2 \left( a -\frac{b}{r^2} \right)' + r \left( a -\frac{b}{r^2} \right) - a r - \frac{b}{r} \\ &=r^2 \frac{2 b}{r^3} + \not{{a r}} - \frac{b}{r} - \not{{a r}} - \frac{b}{r} \\ &=\frac{2 b}{r} - \frac{2 b}{r} \\ &= 0.\end{aligned}

Given the velocity, we can now determine the pressure up to a constant

\begin{aligned}p' &= \frac{\rho}{r} \left( a r + \frac{b}{r} \right)^2 \\ &= \frac{\rho}{r} \left( a^2 r^2 + \frac{b^2}{r^2} + 2 a b \right) \\ &= \rho \left( a^2 r + \frac{b^2}{r^3} + 2 \frac{a b}{r} \right)\end{aligned}


\begin{aligned}p_r -p_0= \rho \left( \frac{1}{{2}} a^2 r^2 - \frac{b^2}{2 r^2} + 2 a b \ln r \right)\end{aligned} \hspace{\stretch{1}}(3.7)

Solution: Part 2. Fixing the constants.

To determine our integration constants we recall that velocity associated with a radial position \mathbf{x} = r \hat{\mathbf{r}} in cylindrical coordinates takes the form

\begin{aligned}\frac{\mathbf{x}}{dt} = \dot{r} \hat{\mathbf{r}} + r \hat{\boldsymbol{\phi}} \dot{\phi},\end{aligned} \hspace{\stretch{1}}(4.8)

where \dot{\phi} is the angular velocity. The cylinder walls therefore have the velocity

\begin{aligned}v = r \dot{\phi},\end{aligned} \hspace{\stretch{1}}(4.9)

so our boundary conditions (given a no-slip assumption for the fluids) are

\begin{aligned}v(R_1) &= R_1 \Omega_1 \\ v(R_2) &= R_2 \Omega_2.\end{aligned} \hspace{\stretch{1}}(4.10)

This gives us a pair of equations to solve for a and b

\begin{aligned}R_1 \Omega_1 &= a R_1 + \frac{b}{R_1} \\ R_2 \Omega_2 &= a R_2 + \frac{b}{R_2}.\end{aligned} \hspace{\stretch{1}}(4.12)

Multipling each by R_1 and R_2 respectively gives us

\begin{aligned}b = R_1^2 (\Omega_1 - a) = R_2^2 (\Omega_2 - a).\end{aligned} \hspace{\stretch{1}}(4.14)

Rearranging for a we find

\begin{aligned}R_1^2 \Omega_1 - R_2^2 \Omega_2 = (R_1^2 - R_2^2) a,\end{aligned} \hspace{\stretch{1}}(4.15)


\begin{aligned}a = \frac{ R_2^2 \Omega_2 - R_1^2 \Omega_1}{R_2^2 - R_1^2}.\end{aligned} \hspace{\stretch{1}}(4.16)

For b we have

\begin{aligned}b &= R_1^2 (\Omega_1 - a) \\ &=\frac{R_1^2 }{R_2^2 - R_1^2}(\Omega_1 ( R_2^2 - \not{{R_1^2}}) - R_2^2 \Omega_2 + \not{{R_1^2 \Omega_1}}),\end{aligned}


\begin{aligned}b = \frac{R_1^2 R_2^2}{R_2^2 - R_1^2} (\Omega_1 -\Omega_2).\end{aligned} \hspace{\stretch{1}}(4.17)

This gives us


\begin{aligned}v(r) = \frac{1}{{R_2^2 - R_1^2}}\left(\left( R_2^2 \Omega_2 - R_1^2 \Omega_1\right) r+\frac{R_1^2 R_2^2}{r} (\Omega_1 -\Omega_2)\right)\end{aligned} \hspace{\stretch{1}}(4.18a)

\begin{aligned}\begin{aligned}p_r -&p_0= \frac{\rho }{(R_2^2 - R_1^2)^2} \times \\ &\left( \frac{1}{{2}} \left( R_2^2 \Omega_2 - R_1^2 \Omega_1\right)^2r^2 -\frac{R_1^4 R_2^4}{2 r^2} (\Omega_1 - \Omega_2)^2+ 2 \left( R_2^2 \Omega_2 - R_1^2 \Omega_1\right) R_1^4 R_2^4 (\Omega_1 - \Omega_2)^2 \ln r\right).\end{aligned}\end{aligned} \hspace{\stretch{1}}(4.18b)


Solution: Part 3. Friction and torque.

We can expand out the identity for the traction vector

\begin{aligned}\mathbf{t} = \mathbf{e}_i \sigma_{ij} n_j= -p \hat{\mathbf{n}} + \mu \left( 2 (\hat{\mathbf{n}} \cdot \boldsymbol{\nabla}) \mathbf{u} + \hat{\mathbf{n}} \times (\boldsymbol{\nabla} \times \mathbf{u})\right),\end{aligned} \hspace{\stretch{1}}(5.19)

in cylindrical coordinates and find


\begin{aligned}\sigma_{rr}=-p + 2 \mu \not{{\frac{\partial {u_r}}{\partial {r}}}}\end{aligned} \hspace{\stretch{1}}(5.20a)

\begin{aligned}\sigma_{\phi \phi}=-p + 2 \mu \left(\frac{1}{{r}}\not{{\frac{\partial {u_\phi}}{\partial {\phi}}}} + \not{{\frac{u_r}{r}}}\right)\end{aligned} \hspace{\stretch{1}}(5.20b)

\begin{aligned}\sigma_{z z}=-p + 2 \mu \not{{\frac{\partial {u_z}}{\partial {z}}}}\end{aligned} \hspace{\stretch{1}}(5.20c)

\begin{aligned}\sigma_{r \phi}=\mu \left( \frac{\partial {u_\phi}}{\partial {r}}+\frac{1}{{r}} \not{{\frac{\partial {u_r}}{\partial {\phi}}}}- \frac{u_\phi}{r}\right)\end{aligned} \hspace{\stretch{1}}(5.20d)

\begin{aligned}\sigma_{\phi z}=\mu \left(\frac{1}{r} \not{{\frac{\partial {u_z}}{\partial {\phi}}}} + \not{{\frac{\partial {u_\phi}}{\partial {z}}}}\right)\end{aligned} \hspace{\stretch{1}}(5.20e)

\begin{aligned}\sigma_{z r}=\mu \left(\not{{\frac{\partial {u_r}}{\partial {z}}}}+ \not{{\frac{\partial {u_z}}{\partial {r}}}}\right),\end{aligned} \hspace{\stretch{1}}(5.20f)


so we have


\begin{aligned}\sigma_{rr} = \sigma_{\phi \phi} = \sigma_{z z} = -p \end{aligned} \hspace{\stretch{1}}(5.21a)

\begin{aligned}\sigma_{\phi z} = \sigma_{z r} = 0\end{aligned} \hspace{\stretch{1}}(5.21b)

\begin{aligned}\sigma_{r \phi} = \mu \left( \frac{\partial {u_\phi}}{\partial {r}} - \frac{u_\phi}{r} \right)\end{aligned} \hspace{\stretch{1}}(5.21c)


We want to expand the last of these

\begin{aligned}\sigma_{r \phi} &= \mu \left( \frac{\partial {u_\phi}}{\partial {r}} - \frac{u_\phi}{r} \right) \\ &= \mu \left( a r + \frac{b}{r}\right)' \\ &= \mu \left( a - \frac{b}{r^2}\right).\end{aligned}

So the traction vector, our force per unit area on the fluid at the inner surface (where the normal is \hat{\mathbf{r}}), is

\begin{aligned}\mathbf{t}_1 = -p \hat{\mathbf{r}} + \mu \left( a - \frac{b}{r^2} \right) \hat{\boldsymbol{\phi}},\end{aligned} \hspace{\stretch{1}}(5.22)

and our torque per unit area from the inner cylinder on the fluid is thus

\begin{aligned}\boldsymbol{\tau}_1 = r \hat{\mathbf{r}} \times \mathbf{t} = r \mu \left( a - \frac{b}{r^2} \right) \hat{\mathbf{z}}.\end{aligned} \hspace{\stretch{1}}(5.23)

Observing that our stress tensors flip sign for an inwards normal, our torque per unit area from the outer cylinder is

\begin{aligned}\boldsymbol{\tau}_2 = r \hat{\mathbf{r}} \times (-\mathbf{t}_1) = -r \mu \left( a - \frac{b}{r^2} \right) \hat{\mathbf{z}}.\end{aligned} \hspace{\stretch{1}}(5.24)

For the complete torque on the fluid due to a strip of width \Delta z the magnitudes of the total torque from each cylinder are respectively

\begin{aligned}\tau_1 = 2 \pi r^2 \Delta z \mu \left( a - \frac{b}{r^2} \right)\end{aligned} \hspace{\stretch{1}}(5.25)

\begin{aligned}\tau_2 = - 2 \pi r^2 \Delta z \mu \left( a - \frac{b}{r^2} \right)\end{aligned} \hspace{\stretch{1}}(5.26)

As expected these torques on the fluids sum to zero

\begin{aligned}\tau_2 + \tau_1 = 0.\end{aligned} \hspace{\stretch{1}}(5.27)

Evaluating these at R_1 and R_2 respectively gives us the torques on the fluid by the cylinders, so inverting these provides the torques on the cylinders by the fluid

\begin{aligned}\text{torque on inner cylinder (1) by the fluid} = -(2 \pi R_1^2) \Delta z \mu \left( a - \frac{b}{R_1^2} \right).\end{aligned} \hspace{\stretch{1}}(5.28)

For the outer cylinder the total torque on a strip of width \Delta z is

\begin{aligned}\text{torque on outer cylinder (2) by the fluid} = (2 \pi R_2^2) \Delta z \mu \left( a - \frac{b}{R_2^2} \right).\end{aligned} \hspace{\stretch{1}}(5.29)


[1] D.J. Acheson. Elementary fluid dynamics. Oxford University Press, USA, 1990.

[2] Wikipedia. Taylor-couette flow — wikipedia, the free encyclopedia [online]. 2012. [Online; accessed 10-April-2012].\%E2\%80\%93Couette_flow&oldid=483281707.

Posted in Math and Physics Learning. | Tagged: , , , , , | Leave a Comment »

Potential due to cylindrical distribution.

Posted by peeterjoot on February 27, 2012

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


Consider a cylindrical distribution of mass (or charge) as in figure (\ref{fig:cylinderPotential:cylinderPotentialFig1}), with points in the cylinder given by \mathbf{r}' = (r', \theta ', z') coordinates, and the point of measurement of the potential measured at \mathbf{r} = (r, 0, 0).

\caption{coordinates for evaluation of cylindrical potential.}

Our potential, for a uniform distribution, will be proportional to

\begin{aligned}\begin{aligned}\phi(r) &= \int \frac{ dV'}{{\left\lvert{\mathbf{r} - \mathbf{r}'}\right\rvert}} \\ &=\int_0^R r' dr' \int_0^{2 \pi } d\theta' \int_{-L}^L \frac{dz' }{\sqrt{(z')^2+ {\left\lvert{r - r' e^{i \theta }}\right\rvert}^2}}.\end{aligned}\end{aligned} \hspace{\stretch{1}}(1.1)

Attempting to evaluate the integrals.


\begin{aligned}\int_{-L}^L \frac{1}{\sqrt{z^2+u^2}} \, dz=\log \left(\frac{L+\sqrt{L^2+u^2}}{-L+\sqrt{L^2+u^2}}\right)\end{aligned} \hspace{\stretch{1}}(2.2)

This is found to be

\begin{aligned}\phi(\mathbf{r}) = \int_0^R r' dr' \int_0^{2 \pi } d\theta'\log\left(\frac{L+\sqrt{L^2+ {\left\lvert{r - r' e^{ i \theta }}\right\rvert}^2}}{-L+\sqrt{L^2+ {\left\lvert{r - r' e^{ i \theta }}\right\rvert}^2}}\right)\end{aligned} \hspace{\stretch{1}}(2.3)

It is clear that we can’t evaluate this limit directly for L \rightarrow \infty since that gives us \infty/0 in the logarithm term. Presuming this can be evaluated, we must have to evaluate the complete set of integrals first, then take the limit. Based on the paper, it appears that this can be evaluated, however, the approach used therein uses mathematics a great deal more sophisticated than I can grasp without a lot of study.

Can we proceed blindly using computational tools to do the work? Attempting to evaluate the remaining integrals with Mathematica fails, since evaluation of both

\begin{aligned}\int_0^R r dr \log \left(a+\sqrt{a^2+r^2+b^2-2 r b \cos (\theta )}\right) \end{aligned} \hspace{\stretch{1}}(2.4)


\begin{aligned}\int_0^{2 \pi } d\theta\log \left(a+\sqrt{a^2+r^2+b^2-2 r b \cos (\theta )}\right) ,\end{aligned} \hspace{\stretch{1}}(2.5)

either time out, or take long enough that I aborted the attempt to let them evaluate.

Alternate evaluation order?

We can also attempt to evaluate this by integrating in different orders. We can for example do the r' coordinate integral first

\begin{aligned}\begin{aligned}\int_0^R &dr\frac{r}{\sqrt{a^2+r^2-2 r b \cos (\theta )+b^2}}  \\ &=\sqrt{a^2+b^2-2 b R \cos (\theta )+R^2} -\sqrt{a^2+b^2} \\ &+b \cos (\theta ) \log \left(\frac{\sqrt{a^2+b^2-2 b R \cos (\theta )+R^2}-b \cos (\theta )+R}{\sqrt{a^2+b^2}-b \cos (\theta )}\right).\end{aligned}\end{aligned} \hspace{\stretch{1}}(2.6)

It should be noted that this returns a number of hard to comprehend ConditionalExpression terms, so care manipulating this expression may also be required.

If we try the angular integral first, we get

\begin{aligned}\int_0^{2 \pi } d\theta\frac{1}{\sqrt{a^2+r^2-2 r b \cos (\theta )+b^2}} =\frac{2 K\left(-\frac{4 b r}{a^2+(b-r)^2}\right)}{\sqrt{a^2+(b-r)^2}}+\frac{2 K\left(\frac{4 b r}{a^2+(b+r)^2}\right)}{\sqrt{a^2+(b+r)^2}}\end{aligned} \hspace{\stretch{1}}(2.7)


\begin{aligned}K(m) = F\left( \left.\frac{\pi }{2}\right| m \right)\end{aligned} \hspace{\stretch{1}}(2.8)

is the complete elliptic integral of the first kind. Actually evaluating this integral, especially in the limiting case, probably requires stepping back and thinking a bit (or a lot) instead of blindly trying to evaluate.

Posted in Math and Physics Learning. | Tagged: , , | Leave a Comment »

Potential for an infinitesimal width infinite plane. Take III

Posted by peeterjoot on February 24, 2012

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

Document generation experiment.

This little document was generated as an experiment using Mathematica and some post processing in latex.

The File menu save as latex produced latex that couldn’t be compiled, but mouse selected, copy-as latex worked out fairly well.

Post processing done included:

\item Adding in latex prologue.
\item Stripping out the text boxes.
\item Adding in equation environments.
\item Latex generation for math output in inline text sections was uniformly poor.


I’d like to attempt again to evaluate the potential for infinite plane distribution. The general form of our potential takes the form

\begin{aligned}\phi(\mathbf{x}) = G \rho \int \frac{1}{{{\left\lvert{\mathbf{x} - \mathbf{x}'}\right\rvert}}} dV'\end{aligned} \hspace{\stretch{1}}(2.1)

We want to evaluate this with cylindrical coordinates (r', \theta', z'), for a width \epsilon, and radius r, at distance z from the plane.

\begin{aligned}\phi (z, \epsilon , r)= 2 \pi  G \sigma  \frac{1}{\epsilon }\int _{r' = 0}^r\int _{z' = 0}^{\epsilon }\frac{r'}{\sqrt{\left(z-z'\right)^2+\left(r'\right)^2}}dz'dr'\end{aligned} \hspace{\stretch{1}}(2.2)

With the assumption that we will take the limits \epsilon \rightarrow 0, and r \rightarrow \infty. With r^2 = c/\epsilon, this does not converge. How about with r = c/\epsilon?

Performing the r’ integration (with r^2 = c/\epsilon) we find

\begin{aligned}\phi (z, \epsilon )= 2 \pi  G \sigma  \frac{1}{\epsilon }\int_{z' = 0}^{\epsilon } \left(\sqrt{\frac{c^2}{\epsilon ^2}+(z-z')^2}-\sqrt{(z-z')^2}\right) \, dz'\end{aligned} \hspace{\stretch{1}}(2.3)

Attempting to let \textit{Mathematica} evaluate this takes a long time. Long enough that I aborted the attempt to evaluate it.

Instead, first evaluating the z’ integral we have

\begin{aligned}\phi (z, \epsilon , r)=\frac{2 \pi  G \sigma }{\epsilon }\int _{r' = 0}^{c/\epsilon }\left(\log \left(\sqrt{\left(r'\right)^2+z^2}+z\right)-\log \left(\sqrt{\left(r'\right)^2+(z-\epsilon )^2}+z-\epsilon \right)\right)dr'\end{aligned} \hspace{\stretch{1}}(2.4)

This second integral can then be evaluated in reasonable time:

\begin{aligned}\begin{aligned}\phi (z, \epsilon )&= \frac{2 \pi  G \sigma }{\epsilon ^2} \left(c \log \left(\frac{\sqrt{\frac{c^2}{\epsilon ^2}+z^2}+z}{\sqrt{\frac{c^2}{\epsilon ^2}+(z-\epsilon )^2}+z-\epsilon }\right)+\epsilon  \left(z \log \left(\frac{(z-\epsilon ) \left(\sqrt{c^2+z^2 \epsilon ^2}+c\right)}{z}\right)+(\epsilon -z) \log \left(\sqrt{c^2+\epsilon ^2 (z-\epsilon )^2}+c\right)-\epsilon  \log (\epsilon  (z-\epsilon ))\right)\right) \\ &=2 \pi  G \sigma  \left(\frac{c}{\epsilon ^2}\log \left(\frac{\sqrt{c^2+z^2 \epsilon ^2}+z \epsilon }{\sqrt{c^2+\epsilon ^2(z-\epsilon )^2}+\epsilon (z-\epsilon )}\right)+\frac{z}{\epsilon } \log \left(\frac{(z-\epsilon ) \left(\sqrt{c^2+z^2 \epsilon ^2}+c\right)}{z\left(\sqrt{c^2+\epsilon ^2 (z-\epsilon )^2}+c\right)}\right)+ \log \left(\frac{\sqrt{c^2+\epsilon ^2 (z-\epsilon )^2}+c}{\epsilon  (z-\epsilon )}\right)\right)\end{aligned}\end{aligned} \hspace{\stretch{1}}(2.5)

Does this have a limit as \epsilon \rightarrow 0? No, the last term is clearly divergent for c \neq 0.

Posted in Math and Physics Learning. | Tagged: , , , , | Leave a Comment »

PHY454H1S Continuum Mechanics. Lecture 2. Introduction and strain tensor. Taught by Prof. K. Das.

Posted by peeterjoot on January 21, 2012

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


Peeter’s lecture notes from class. May not be entirely coherent.


Mechanics could be defined as the study of effects of forces and displacements on a physical body

\caption{Physical body.}

In continuum mechanics we have a physical body and we are interested in the internal motions in the object.

\caption{Control volume elements.}

For the first time considering mechanics we have to introduce the concepts of fields to make progress tackling these problems.

We will have use of the following types of fields

\item Scalar fields. 3^0 components. Examples: density, Temperature, …
\item Vector fields. 3^1 components. Examples: Force, velocity.
\item Tensor fields. 3^2 components. Examples: stress, strain.

We have to consider objects (a control volume) that is small enough that we can consider that we have a point in space limit for the quantities of density and velocity. At the same time we cannot take this limiting process to the extreme, since if we use a control volume that is sufficiently small, quantum and inter-atomic effects would have to be considered.

\caption{Mass and volume ratios at different scales.}

Stress and Strain definitions.


Measure of the Internal force on the surfaces.


Measure of the deformation of the body.

Strain Tensor.

This follows [1] section 1 very closely.

Utilizing summation convention consider a set of small internal displacements u_1, u_2, u_3 to the x, y, z coordinates so that the transformation x_i \rightarrow x_i' is related by

\caption{Differential change to the object.}

\begin{aligned}u_i &= x_i' - x_i \\ x_i' &= g(x_i) \\ u_i &= f(x_i)\end{aligned} \hspace{\stretch{1}}(3.1)

(ie: x_i' is a function of all the initial coordinates, as are the displacements u_i).

\begin{aligned}dx_i' = dx_i + du_i\end{aligned} \hspace{\stretch{1}}(3.4)

\begin{aligned}dl &= \sqrt{dx_k dx_k} \\ dl' &= \sqrt{d{x'}_k d{x'}_k}\end{aligned} \hspace{\stretch{1}}(3.5)


\begin{aligned}{dl'}^2 = (dx_k + du_k)(dx_k + du_k)= dl^2 + 2 dx_k du_k + du_k du_k\end{aligned} \hspace{\stretch{1}}(3.7)


\begin{aligned}du_i = \frac{\partial {u_i}}{\partial {x_k}} dx_k\end{aligned} \hspace{\stretch{1}}(3.8)

we have

\begin{aligned}du_i^2 = \frac{\partial {u_i}}{\partial {x_k}} dx_k\frac{\partial {u_i}}{\partial {x_l}} dx_l\end{aligned} \hspace{\stretch{1}}(3.9)

\begin{aligned}{dl'}^2 &= dl^2 + 2 \frac{\partial {u_i}}{\partial {x_k}} dx_k dx_i + \frac{\partial {u_l}}{\partial {x_i}} \frac{\partial {u_l}}{\partial {x_k}} dx_i dx_k \\ &= dl^2 + \left(\frac{\partial {u_i}}{\partial {x_k}} +\frac{\partial {u_k}}{\partial {x_i}} \right)dx_k dx_i + \frac{\partial {u_l}}{\partial {x_i}} \frac{\partial {u_l}}{\partial {x_k}} dx_i dx_k \\ &=dl^2 + 2 e_{ik} dx_i dx_k\end{aligned}

We write

\begin{aligned}{dl'}^2 - dl^2 = 2 e_{ik} dx_i dx_k\end{aligned} \hspace{\stretch{1}}(3.10)

where we define the \emph{strain tensor} as

\begin{aligned}e_{ik} = \frac{1}{{2}} \left(\left(\frac{\partial {u_i}}{\partial {x_k}} +\frac{\partial {u_k}}{\partial {x_i}} \right)+ \frac{\partial {u_l}}{\partial {x_i}} \frac{\partial {u_l}}{\partial {x_k}} \right)\end{aligned} \hspace{\stretch{1}}(3.11)

Here e_{ik} is a 3 \times 3 matrix in Cartesian coordinates

\begin{aligned}\begin{bmatrix}e_{11} & e_{12} & e_{13} \\ e_{21} & e_{22} & e_{23} \\ e_{31} & e_{32} & e_{33} \\ \end{bmatrix}\end{aligned} \hspace{\stretch{1}}(3.12)

We see from 3.11 that e_{ik} is symmetric, so we have

\begin{aligned}e_{21} &= e_{12} \\ e_{31} &= e_{13} \\ e_{32} &= e_{23}\end{aligned} \hspace{\stretch{1}}(3.13)

Because any real symmetric matrix can be diagonalized we can write in some coordinate system

\begin{aligned}\bar{e}_{ik} = \begin{bmatrix}\bar{e}_{11} & 0 & 0 \\ 0 & \bar{e}_{22} & 0 \\ 0 & 0 & \bar{e}_{33} \\ \end{bmatrix}\end{aligned} \hspace{\stretch{1}}(3.16)

\begin{aligned}{dx_1'}^2 &= (1 + 2 \bar{e}_{11}) dx_1^2 \\ {dx_2'}^2 &= (1 + 2 \bar{e}_{22}) dx_2^2 \\ {dx_3'}^2 &= (1 + 2 \bar{e}_{33}) dx_3^2\end{aligned} \hspace{\stretch{1}}(3.17)

If our changes are small enough we can also write approximately, taking the first order term in the square root evaluation

\begin{aligned}dx_1' &\approx (1 + \bar{e}_{11}) dx_1 \\ dx_2' &\approx (1 + \bar{e}_{22}) dx_2 \\ dx_3' &\approx (1 + \bar{e}_{33}) dx_3\end{aligned} \hspace{\stretch{1}}(3.20)

We are also free to define a volume element

\begin{aligned}dV' = dx_1'dx_2'dx_3'\approx(1 + e_{11})(1 + e_{22})(1 + e_{33})dx_1 dx_2 dx_3\end{aligned} \hspace{\stretch{1}}(3.23)

\begin{aligned}dV' = (1 + e_{11} +e_{22} +e_{33} ) dV\end{aligned} \hspace{\stretch{1}}(3.24)

So the change of volume is given by the trace

\begin{aligned}dV' = ( 1 + e_{ii} )^2 dV\end{aligned} \hspace{\stretch{1}}(3.25)

Strain Tensor in cylindrical coordinates.

At the end of the section in the text, the formulas for the spherical and cylindrical versions (to first order) of the strain tensor is given without derivation. Let’s do that derivation for the cylindrical case, which is simpler. It appears that use of explicit vector notation is helpful here, so we write

\begin{aligned}\mathbf{x} &= r \hat{\mathbf{r}} + z \hat{\mathbf{z}} \\ \mathbf{u} & u_r \hat{\mathbf{r}} + u_\phi \hat{\boldsymbol{\phi}} + u_z \hat{\mathbf{z}}\end{aligned} \hspace{\stretch{1}}(3.26)


\begin{aligned}\hat{\mathbf{r}} &= \mathbf{e}_1 e^{i\phi} \\ \hat{\boldsymbol{\phi}} &= \mathbf{e}_2 e^{i\phi} \\ i &= \mathbf{e}_1 \mathbf{e}_2\end{aligned} \hspace{\stretch{1}}(3.28)

Since \hat{\mathbf{r}} and \hat{\boldsymbol{\phi}} are functions of position, we will need their differentials

\begin{aligned}d\hat{\mathbf{r}} &= \mathbf{e}_1 \mathbf{e}_1 \mathbf{e}_2 e^{i\phi} d\phi = \mathbf{e}_2 e^{i \phi} d\phi \\ d\hat{\boldsymbol{\phi}} &= \mathbf{e}_2 \mathbf{e}_1 \mathbf{e}_2 e^{i\phi} d\phi = -\mathbf{e}_2 e^{i \phi} d\phi,\end{aligned} \hspace{\stretch{1}}(3.31)

but these are just scaled basis vectors

\begin{aligned}d\hat{\mathbf{r}} &= \hat{\boldsymbol{\phi}} d\phi \\ d\hat{\boldsymbol{\phi}} &= -\hat{\mathbf{r}} d\phi.\end{aligned} \hspace{\stretch{1}}(3.33)

So for our \mathbf{x} and \mathbf{u} differentials we find

\begin{aligned}d\mathbf{x} &= dr \hat{\mathbf{r}} + r d\hat{\mathbf{r}} + dz \hat{\mathbf{z}} \\ &= dr \hat{\mathbf{r}} + r \hat{\boldsymbol{\phi}} d\phi + dz \hat{\mathbf{z}},\end{aligned}


\begin{aligned}d\mathbf{u} &= du_r \hat{\mathbf{r}} + du_\phi \hat{\boldsymbol{\phi}} + du_z \hat{\mathbf{z}} + u_r \hat{\boldsymbol{\phi}} d\phi - u_\phi \hat{\mathbf{r}} d\phi \\ &= \hat{\mathbf{r}}( du_r - u_\phi d\phi )+ \hat{\boldsymbol{\phi}} ( du_\phi + u_r d\phi )+ \hat{\mathbf{z}} ( du_z ).\end{aligned}

Putting these together we have

\begin{aligned}d\mathbf{l}' &= d\mathbf{u} + d\mathbf{x} \\ &= \hat{\mathbf{r}}( du_r - u_\phi d\phi + dr )+ \hat{\boldsymbol{\phi}} ( du_\phi + u_r d\phi + r d\phi )+ \hat{\mathbf{z}} ( du_z + dz ).\end{aligned}

For the squared magnitude’s difference from d\mathbf{x}^2 we have

\begin{aligned}(d\mathbf{l}')^2 - d\mathbf{x}^2&= ( du_r - u_\phi d\phi + dr )^2+ ( du_\phi + u_r d\phi + r d\phi )^2+ ( du_z + dz )^2-dr^2 - r^2 d\phi^2 - dz^2 \\ &=( du_r - u_\phi d\phi )^2 + 2 dr ( du_r - u_\phi d\phi )+ ( du_\phi + u_r d\phi )^2+ 2 r d\phi ( du_\phi + u_r d\phi )+ du_z^2 + 2 du_z dz \\ \end{aligned}

Expanding this out, but dropping all the terms that are quadratic in the components of \mathbf{u} or its differentials, we have

\begin{aligned}(d\mathbf{l}')^2 - d\mathbf{x}^2&\approx  2 dr ( du_r - u_\phi d\phi )+ 2 r d\phi ( du_\phi + u_r d\phi )+ 2 du_z dz \\ &=  2 dr du_r - 2 dr u_\phi d\phi + 2 r d\phi du_\phi + 2 r d\phi u_r d\phi + 2 du_z dz \\ &=  2 dr \left( \frac{\partial {u_r}}{\partial {r}} dr+\frac{\partial {u_r}}{\partial {\phi}} d\phi+\frac{\partial {u_r}}{\partial {z}} dz\right) \\ &- 2 dr d\phi u_\phi  \\ &+ 2 r d\phi \left( \frac{\partial {u_\phi}}{\partial {r}} dr+\frac{\partial {u_\phi}}{\partial {\phi}} d\phi+\frac{\partial {u_\phi}}{\partial {z}} dz\right) \\ &+ 2 r d\phi d\phi u_r \\ &+ 2 dz \left( \frac{\partial {u_z}}{\partial {r}} dr+\frac{\partial {u_z}}{\partial {\phi}} d\phi+\frac{\partial {u_z}}{\partial {z}} dz\right) \\ \end{aligned}

Grouping all terms, with all the second order terms neglected, we have

\begin{aligned}\begin{aligned}(d\mathbf{l}')^2 - d\mathbf{x}^2&=2 dr dr \frac{\partial {u_r}}{\partial {r}} + 2 r^2 d\phi d\phi \left( \frac{1}{{r}} \frac{\partial {u_\phi}}{\partial {\phi}} +\frac{1}{{r}} u_r \right)+ 2 dz dz \frac{\partial {u_z}}{\partial {z}}  \\ &+ 2 dz dr \left( \frac{\partial {u_r}}{\partial {z}} + \frac{\partial {u_z}}{\partial {r}} \right)+ 2 dr r d\phi \left( \frac{\partial {u_\phi}}{\partial {r}} - \frac{1}{{r}} u_\phi + \frac{1}{{r}} \frac{\partial {u_r}}{\partial {\phi}} \right)+ 2 dz r d\phi \left( \frac{\partial {u_\phi}}{\partial {z}} +\frac{1}{{r}} \frac{\partial {u_z}}{\partial {\phi}} \right).\end{aligned}\end{aligned} \hspace{\stretch{1}}(3.35)

From this we can read off the result quoted in the text

\begin{aligned}2 e_{rr} &= \frac{\partial {u_r}}{\partial {r}}  \\ 2 e_{\phi\phi} &= \frac{1}{{r}} \frac{\partial {u_\phi}}{\partial {\phi}} +\frac{1}{{r}} u_r  \\ 2 e_{zz} &= \frac{\partial {u_z}}{\partial {z}}  \\ 2 e_{zr} &= \frac{\partial {u_r}}{\partial {z}} + \frac{\partial {u_z}}{\partial {r}} \\ 2 e_{r\phi} &= \frac{\partial {u_\phi}}{\partial {r}} - \frac{1}{{r}} u_\phi + \frac{1}{{r}} \frac{\partial {u_r}}{\partial {\phi}} \\ 2 e_{\phi z} &= \frac{\partial {u_\phi}}{\partial {z}} +\frac{1}{{r}} \frac{\partial {u_z}}{\partial {\phi}}.\end{aligned} \hspace{\stretch{1}}(3.36)

Observe that we have to introduce factors of r along with all the d\phi‘s, when we factored out the tensor components. That’s an important looking detail, which isn’t obvious unless one works through the derivation.

Note that in class we retained the second order terms. That becomes a messier calculation and I’ve cheated using the symbolic capabilities of mathematica to do it

\begin{aligned}\begin{aligned}&(d\mathbf{l}')^2 - d\mathbf{x}^2 \\ &= (dr)^2 \left(2 \frac{\partial u_r}{\partial r}+\left(\frac{\partial u_r}{\partial r}\right)^2+\left(\frac{\partial u_z}{\partial r}\right)^2+\left(\frac{\partial u_{\phi }}{\partial r}\right)^2\right) \\ &+(d\phi )^2 \left(2 r u_r+u_r^2+u_{\phi }^2-2 u_{\phi } \frac{\partial u_r}{\partial \phi }+\left(\frac{\partial u_r}{\partial \phi }\right)^2+\left(\frac{\partial u_z}{\partial \phi }\right)^2+2 r \frac{\partial u_{\phi }}{\partial \phi }+2 u_r \frac{\partial u_{\phi }}{\partial \phi }+\left(\frac{\partial u_{\phi }}{\partial \phi }\right)^2\right) \\ &+(dz)^2 \left(\left(\frac{\partial u_r}{\partial z}\right)^2+2 \frac{\partial u_z}{\partial z}+\left(\frac{\partial u_z}{\partial z}\right)^2+\left(\frac{\partial u_{\phi }}{\partial z}\right)^2\right) \\ &+dr d\phi  \left(-2 u_{\phi }-2 u_{\phi } \frac{\partial u_r}{\partial r}+2 \frac{\partial u_r}{\partial \phi }+2 \frac{\partial u_r}{\partial r} \frac{\partial u_r}{\partial \phi }+2 \frac{\partial u_z}{\partial r} \frac{\partial u_z}{\partial \phi }+2 r \frac{\partial u_{\phi }}{\partial r}+2 u_r \frac{\partial u_{\phi }}{\partial r}+2 \frac{\partial u_{\phi }}{\partial r} \frac{\partial u_{\phi }}{\partial \phi }\right) \\ &+dz d\phi  \left(-2 u_{\phi } \frac{\partial u_r}{\partial z}+2 \frac{\partial u_r}{\partial z} \frac{\partial u_r}{\partial \phi }+2 \frac{\partial u_z}{\partial \phi }+2 \frac{\partial u_z}{\partial z} \frac{\partial u_z}{\partial \phi }+2 r \frac{\partial u_{\phi }}{\partial z}+2 u_r \frac{\partial u_{\phi }}{\partial z}+2 \frac{\partial u_{\phi }}{\partial z} \frac{\partial u_{\phi }}{\partial \phi }\right) \\ &+dr dz \left(2 \frac{\partial u_r}{\partial z}+2 \frac{\partial u_r}{\partial r} \frac{\partial u_r}{\partial z}+2 \frac{\partial u_z}{\partial r}+2 \frac{\partial u_z}{\partial r} \frac{\partial u_z}{\partial z}+2 \frac{\partial u_{\phi }}{\partial r} \frac{\partial u_{\phi }}{\partial z}\right).\end{aligned}\end{aligned} \hspace{\stretch{1}}(3.42)

As with the first order case, we can read off the tensor coordinates by inspection (once we factor out the various factors of 2 and r). The next logical step would be to do the spherical tensor calculation. That would likely be particularily messy if we attempted it in the brute force fashion. Let’s step back and look at the general case, before tackling there sphereical polar form explicitly.

Strain Tensor for general coordinate representation.

Now let’s dispense with the assumption that we have an orthonormal frame. Given an arbitrary, not neccessarily orthonormal, position dependent frame \{e_\mu\}, and its reciprocal frame \{e^\mu\}, as defined by

\begin{aligned}e_\mu \cdot e^\nu = {\delta_\mu}^\nu.\end{aligned} \hspace{\stretch{1}}(3.43)

Our coordinate representation, with summation and dimensionality implied, is

\begin{aligned}\mathbf{x} &= x^\mu e_\mu = x_\nu e^\nu \\ \mathbf{u} &= u^\mu e_\mu = u_\nu e^\nu.\end{aligned} \hspace{\stretch{1}}(3.44)

Our differentials are

\begin{aligned}\begin{aligned}d\mathbf{x} &= dx^\mu e_\mu + x^\mu d e_\mu \\ &= \sum_\alpha d\alpha \left( \frac{\partial {x^\mu}}{\partial {\alpha}} e_\mu+x^\mu\frac{\partial {e_\mu}}{\partial {\alpha}} \right),\end{aligned}\end{aligned} \hspace{\stretch{1}}(3.46)


\begin{aligned}\begin{aligned}d\mathbf{u} &= du^\mu e_\mu + u^\mu d e_\mu \\ &= \sum_\alpha d\alpha \left( \frac{\partial {u^\mu}}{\partial {\alpha}} e_\mu+u^\mu\frac{\partial {e_\mu}}{\partial {\alpha}} \right).\end{aligned}\end{aligned} \hspace{\stretch{1}}(3.47)

Summing these we have

\begin{aligned}d\mathbf{u} + d\mathbf{u} = \sum_\alpha d\alpha \left( \left(\frac{\partial {x^\mu}}{\partial {\alpha}} +\frac{\partial {u^\mu}}{\partial {\alpha}} \right)e_\mu+\left(x^\mu+u^\mu\right)\frac{\partial {e_\mu}}{\partial {\alpha}} \right).\end{aligned} \hspace{\stretch{1}}(3.48)

Taking dot products to form the squares we have

\begin{aligned}d\mathbf{x}^2 &= \sum_{\alpha, \beta} d\alpha d\beta \left( \frac{\partial {x^\mu}}{\partial {\alpha}} e_\mu+x^\mu\frac{\partial {e_\mu}}{\partial {\alpha}} \right)\cdot\left( \frac{\partial {x_\nu}}{\partial {\beta}} e^\nu+x_\nu\frac{\partial {e^\nu}}{\partial {\beta}} \right) \\ &=\sum_{\alpha, \beta} d\alpha d\beta \left( \frac{\partial {x^\mu}}{\partial {\alpha}} \frac{\partial {x_\mu}}{\partial {\beta}} +x^\mu x_\nu\frac{\partial {e_\mu}}{\partial {\alpha}} \cdot\frac{\partial {e^\nu}}{\partial {\beta}} + 2 \frac{\partial {x^\mu}}{\partial {\alpha}} x_\nu e_\mu \cdot\frac{\partial {e^\nu}}{\partial {\beta}} \right),\end{aligned}


\begin{aligned}&(d\mathbf{u} + d\mathbf{x})^2 \\ &= \sum_{\alpha, \beta}d\alpha d\beta \left( \left(\frac{\partial {x^\mu}}{\partial {\alpha}} +\frac{\partial {u^\mu}}{\partial {\alpha}} \right)e_\mu+\left(x^\mu+u^\mu\right)\frac{\partial {e_\mu}}{\partial {\alpha}} \right)\cdot\left( \left(\frac{\partial {x_\nu}}{\partial {\beta}} +\frac{\partial {u_\nu}}{\partial {\beta}} \right)e^\nu+\left(x_\nu+u_\nu\right)\frac{\partial {e^\nu}}{\partial {\beta}} \right) \\ &= \sum_{\alpha, \beta}d\alpha d\beta \left(\left(\frac{\partial {x^\mu}}{\partial {\alpha}} +\frac{\partial {u^\mu}}{\partial {\alpha}} \right)\left(\frac{\partial {x_\mu}}{\partial {\beta}} +\frac{\partial {u_\mu}}{\partial {\beta}} \right)+\left(x^\mu+u^\mu\right)\left(x_\nu+u_\nu\right)\frac{\partial {e_\mu}}{\partial {\alpha}} \cdot\frac{\partial {e^\nu}}{\partial {\beta}} +2\left(x^\mu+u^\mu\right)e^\nu\cdot\frac{\partial {e_\mu}}{\partial {\alpha}} \left(\frac{\partial {x_\nu}}{\partial {\beta}} +\frac{\partial {u_\nu}}{\partial {\beta}} \right)\right).\end{aligned}

Taking the difference we find

\begin{aligned}\begin{aligned}&(d\mathbf{u} + d\mathbf{x})^2 - d\mathbf{x}^2 \\ &=\sum_{\alpha, \beta}d\alpha d\beta \left( \frac{\partial {u^\mu}}{\partial {\alpha}} \frac{\partial {u_\mu}}{\partial {\beta}} +2\frac{\partial {u^\mu}}{\partial {\alpha}} \frac{\partial {x_\mu}}{\partial {\beta}} + \left(u^\mu u_\nu +x^\mu u_\nu +u^\mu x_\nu \right)\frac{\partial {e_\mu}}{\partial {\alpha}}\cdot\frac{\partial {e^\nu}}{\partial {\beta}} +2 \left(\frac{\partial {x^\mu}}{\partial {\alpha}}u_\nu+\frac{\partial {u^\mu}}{\partial {\alpha}}(x_\nu+u_\nu)\right)e_\mu \cdot \frac{\partial {e^\nu}}{\partial {\beta}}\right).\end{aligned}\end{aligned} \hspace{\stretch{1}}(3.49)

To evaluate this, it is useful, albeit messier, to group terms a bit

\begin{aligned}\begin{aligned}&(d\mathbf{u} + d\mathbf{x})^2 - d\mathbf{x}^2 \\ &=\sum_{\alpha}2 d\alpha d\alpha \left( \frac{1}{{2}}\frac{\partial {u^\mu}}{\partial {\alpha}} \frac{\partial {u_\mu}}{\partial {\alpha}} +\frac{\partial {u^\mu}}{\partial {\alpha}} \frac{\partial {x_\mu}}{\partial {\alpha}} + \frac{1}{{2}}\left(u^\mu u_\nu +x^\mu u_\nu +u^\mu x_\nu \right)\frac{\partial {e_\mu}}{\partial {\alpha}}\cdot\frac{\partial {e^\nu}}{\partial {\alpha}} +\left(\frac{\partial {x^\mu}}{\partial {\alpha}}u_\nu+\frac{\partial {u^\mu}}{\partial {\alpha}}(x_\nu+u_\nu)\right)e_\mu \cdot \frac{\partial {e^\nu}}{\partial {\alpha}}\right) \\ &+\sum_{\alpha < \beta}2 d\alpha d\beta \left( \frac{\partial {u^\mu}}{\partial {\alpha}} \frac{\partial {u_\mu}}{\partial {\beta}} +\frac{\partial {u^\mu}}{\partial {\alpha}} \frac{\partial {x_\mu}}{\partial {\beta}} +\frac{\partial {u^\mu}}{\partial {\beta}} \frac{\partial {x_\mu}}{\partial {\alpha}} + \frac{1}{{2}}\left(u^\mu u_\nu +x^\mu u_\nu +u^\mu x_\nu \right)\left(\frac{\partial {e_\mu}}{\partial {\alpha}}\cdot\frac{\partial {e^\nu}}{\partial {\beta}} +\frac{\partial {e_\mu}}{\partial {\beta}}\cdot\frac{\partial {e^\nu}}{\partial {\alpha}} \right) \right) \\ &+\sum_{\alpha < \beta}2 d\alpha d\beta \left( \left(\frac{\partial {x^\mu}}{\partial {\alpha}}u_\nu+\frac{\partial {u^\mu}}{\partial {\alpha}}(x_\nu+u_\nu)\right)e_\mu \cdot \frac{\partial {e^\nu}}{\partial {\beta}}+\left(\frac{\partial {x^\mu}}{\partial {\beta}}u_\nu+\frac{\partial {u^\mu}}{\partial {\beta}}(x_\nu+u_\nu)\right)e_\mu \cdot \frac{\partial {e^\nu}}{\partial {\alpha}}\right)\end{aligned}\end{aligned} \hspace{\stretch{1}}(3.50)

Here \alpha < \beta is used to denote summation over the pairs \alpha \ne \beta just once, not neccessarily any numeric ordering. For example with \alpha, \beta \in \{r, \phi, z\}, this could be the set \{\alpha, \beta\} \in \{r \phi, \phi z, z r\}.

Cartesian tensor.

In the Cartesian case all the partials of the unit vectors are zero, and we also have no need of upper or lower indexes. We are left with just

\begin{aligned}(d\mathbf{u} + d\mathbf{x})^2 - d\mathbf{x}^2 =\sum_{i, j, k}dx^idx^j\left( \frac{\partial {u^k}}{\partial {x^i}} \frac{\partial {u^k}}{\partial {x^j}} +2\frac{\partial {u^k}}{\partial {x^i}} \frac{\partial {x^k}}{\partial {x^j}} \right)\end{aligned} \hspace{\stretch{1}}(3.51)

However, since we also have {\partial {x^k}}/{\partial {x^j}} = \delta_{jk}, this is

\begin{aligned}(d\mathbf{u} + d\mathbf{x})^2 - d\mathbf{x}^2 =\sum_{i, j}2dx^idx^j\left( \frac{1}{{2}}\sum_k\frac{\partial {u^k}}{\partial {x^i}} \frac{\partial {u^k}}{\partial {x^j}} +\frac{\partial {u^j}}{\partial {x^i}} \right).\end{aligned} \hspace{\stretch{1}}(3.52)

This essentially recovers the result 3.11 derived in class.

Cylindrial tensor.

Now lets do the cylindrical tensor again, but this time without resorting mathematica brute force.

First we recall that all our basis vector derivatives are zero except for the \phi derivatives, and for those we have

\begin{aligned}\frac{\partial {\hat{\mathbf{r}}}}{\partial {\phi}} &= \hat{\boldsymbol{\phi}} \\ \frac{\partial {\hat{\boldsymbol{\theta}}}}{\partial {\phi}} &= -\hat{\mathbf{r}}.\end{aligned} \hspace{\stretch{1}}(3.53)

If we write

\begin{aligned}\mathbf{x} = r \hat{\mathbf{r}} + z \hat{\mathbf{z}} = x_r \hat{\mathbf{r}} + x_\phi \hat{\boldsymbol{\phi}} + x_z \hat{\mathbf{z}}\end{aligned} \hspace{\stretch{1}}(3.55)

We have for all the x^\mu partials

\begin{aligned}\frac{\partial {x^\mu}}{\partial {\alpha}} = \left\{\begin{array}{l l}1 & \quad \mbox{if latex \alpha = x^\mu = r$ or \alpha = x^\mu = z} \\ 0 & \quad \mbox{otherwise}\end{array}\right.\end{aligned} \hspace{\stretch{1}}(3.56)$

We are now set to evaluate the terms in the sum of 3.50 for the cylindrical coordinate system and shouldn’t need Mathematica to do it. Let’s do this one at a time, starting with all the squared differential pairs. Those are, for \alpha \in \{r, \phi, z\} the value of

\begin{aligned}2 d\alpha d\alpha \left( \frac{1}{{2}}\frac{\partial {u_m}}{\partial {\alpha}} \frac{\partial {u_m}}{\partial {\alpha}} +\frac{\partial {u_m}}{\partial {\alpha}} \frac{\partial {x_m}}{\partial {\alpha}} + \frac{1}{{2}}\left(u_m u_n +x_m u_n +u_m x_n \right)\frac{\partial {e_m}}{\partial {\alpha}}\cdot\frac{\partial {e_n}}{\partial {\alpha}} +\left(\frac{\partial {x_m}}{\partial {\alpha}}u_n+\frac{\partial {u_m}}{\partial {\alpha}}(x_n+u_n)\right)e_m \cdot \frac{\partial {e_n}}{\partial {\alpha}}\right)\end{aligned} \hspace{\stretch{1}}(3.60)

For both r and z all our unit vectors have zero derivatives so we are left respectively with

\begin{aligned}2 dr dr \left( \frac{1}{{2}}\frac{\partial {u_m}}{\partial {r}} \frac{\partial {u_m}}{\partial {r}} +\frac{\partial {u_r}}{\partial {r}} \right),\end{aligned} \hspace{\stretch{1}}(3.60)


\begin{aligned}2 dz dz \left( \frac{1}{{2}}\frac{\partial {u_m}}{\partial {z}} \frac{\partial {u_m}}{\partial {z}} +\frac{\partial {u_z}}{\partial {z}} \right).\end{aligned} \hspace{\stretch{1}}(3.60)

For the \alpha = \phi term we have

\begin{aligned}&2 d\phi d\phi \left( \frac{1}{{2}}\frac{\partial {u_m}}{\partial {\phi}} \frac{\partial {u_m}}{\partial {\phi}} + \frac{1}{{2}}\sum_{m = r, \phi}\left(u_m u_m +2 x_m u_m \right)+\sum_{m n \in \{r \phi, \phi r\}}\left(\frac{\partial {x_m}}{\partial {\phi}}u_n+\frac{\partial {u_m}}{\partial {\phi}}(x_n+u_n)\right)e_m \cdot \frac{\partial {e_n}}{\partial {\phi}}\right) \\ &=2 d\phi d\phi \left( \frac{1}{{2}}\frac{\partial {u_m}}{\partial {\phi}} \frac{\partial {u_m}}{\partial {\phi}} + \frac{1}{{2}} \left( u_r^2 + u_\phi^2 \right) + r u_r-\frac{\partial {u_r}}{\partial {\phi}}u_\phi+\frac{\partial {u_\phi}}{\partial {\phi}}(r+u_r)\right)\end{aligned}

Now, on to the mixed terms. The easiest is the dz dr term, for which all the unit vector derivatives are zero, and we are left with just

\begin{aligned}2 dz dr \left( \frac{\partial {u_m}}{\partial {z}} \frac{\partial {u_m}}{\partial {r}} +\frac{\partial {u_m}}{\partial {z}} \frac{\partial {x_m}}{\partial {r}} +\frac{\partial {u_m}}{\partial {r}} \frac{\partial {x_m}}{\partial {z}} \right)=2 dz dr \left( \frac{\partial {u_m}}{\partial {z}} \frac{\partial {u_m}}{\partial {r}} +\frac{\partial {u_r}}{\partial {z}} +\frac{\partial {u_z}}{\partial {r}} \right)\end{aligned}

Now we have the two messy mixed terms. For the r, \phi term we get

\begin{aligned}&2 dr d\phi \left( \frac{\partial {u_m}}{\partial {r}} \frac{\partial {u_m}}{\partial {\phi}} +\frac{\partial {u_m}}{\partial {r}} \not{{\frac{\partial {x_m}}{\partial {\phi}}}}+\frac{\partial {u_m}}{\partial {\phi}} \frac{\partial {x_m}}{\partial {r}} + \frac{1}{{2}}\left(u_m u_n +x_m u_n +u_m x_n \right)\left(\not{{\frac{\partial {e_m}}{\partial {r}}}}\cdot\frac{\partial {e_n}}{\partial {\phi}} +\frac{\partial {e_m}}{\partial {\phi}}\cdot\not{{\frac{\partial {e_n}}{\partial {r}} }}\right) \right) \\ &+2 dr d\phi \left( \left(\frac{\partial {x_m}}{\partial {r}}u_n+\frac{\partial {u_m}}{\partial {r}}(x_n+u_n)\right)e_m \cdot \frac{\partial {e_n}}{\partial {\phi}}+\left(\frac{\partial {x_m}}{\partial {\phi}}u_n+\frac{\partial {u_m}}{\partial {\phi}}(x_n+u_n)\right)e_m \cdot \not{{\frac{\partial {e_n}}{\partial {r}}}}\right) \\ &=2 dr d\phi \left( \frac{\partial {u_m}}{\partial {r}} \frac{\partial {u_m}}{\partial {\phi}} +\frac{\partial {u_r}}{\partial {\phi}} +u_n\hat{\mathbf{r}} \cdot \frac{\partial {e_n}}{\partial {\phi}}+\frac{\partial {u_m}}{\partial {r}}(x_n+u_n)e_m \cdot \frac{\partial {e_n}}{\partial {\phi}}\right) \\ &=2 dr d\phi \left( \frac{\partial {u_m}}{\partial {r}} \frac{\partial {u_m}}{\partial {\phi}} +\frac{\partial {u_r}}{\partial {\phi}} -u_\phi+\frac{\partial {u_r}}{\partial {r}}(x_n+u_n)\hat{\mathbf{r}} \cdot \frac{\partial {e_n}}{\partial {\phi}}+\frac{\partial {u_\phi}}{\partial {r}}(x_n+u_n)\hat{\boldsymbol{\phi}} \cdot \frac{\partial {e_n}}{\partial {\phi}}\right) \\ &=2 dr d\phi \left( \frac{\partial {u_m}}{\partial {r}} \frac{\partial {u_m}}{\partial {\phi}} +\frac{\partial {u_r}}{\partial {\phi}} -u_\phi-\frac{\partial {u_r}}{\partial {r}}u_\phi+\frac{\partial {u_\phi}}{\partial {r}}(r +u_r)\right) \\ \end{aligned}

Finally for the z, \phi term we have

\begin{aligned}&2 dz d\phi \left( \frac{\partial {u_m}}{\partial {z}} \frac{\partial {u_m}}{\partial {\phi}} +\frac{\partial {u_m}}{\partial {z}} \not{{\frac{\partial {x_m}}{\partial {\phi}} }}+\frac{\partial {u_m}}{\partial {\phi}} \frac{\partial {x_m}}{\partial {z}} + \frac{1}{{2}}\left(u_m u_n +x_m u_n +u_m x_n \right)\left(\not{{\frac{\partial {e_m}}{\partial {z}}}}\cdot\frac{\partial {e_n}}{\partial {\phi}} +\frac{\partial {e_m}}{\partial {\phi}}\cdot\not{{\frac{\partial {e_n}}{\partial {z}} }}\right) \right) \\ &+2 d\phi dz \left( \left(\frac{\partial {x_m}}{\partial {z}}u_n+\frac{\partial {u_m}}{\partial {z}}(x_n+u_n)\right)e_m \cdot \frac{\partial {e_n}}{\partial {\phi}}+\left(\frac{\partial {x_m}}{\partial {\phi}}u_n+\frac{\partial {u_m}}{\partial {\phi}}(x_n+u_n)\right)e_m \cdot \not{{\frac{\partial {e_n}}{\partial {z}}}}\right) \\ &=2 dz d\phi \left( \frac{\partial {u_m}}{\partial {z}} \frac{\partial {u_m}}{\partial {\phi}} +\frac{\partial {u_m}}{\partial {\phi}} \frac{\partial {x_m}}{\partial {z}} +\not{{u_n\hat{\mathbf{z}} \cdot \frac{\partial {e_n}}{\partial {\phi}}}}+\frac{\partial {u_m}}{\partial {z}}(x_n+u_n)e_m \cdot \frac{\partial {e_n}}{\partial {\phi}}\right) \\ &=2 dz d\phi \left( \frac{\partial {u_m}}{\partial {z}} \frac{\partial {u_m}}{\partial {\phi}} +\frac{\partial {u_z}}{\partial {\phi}} -\frac{\partial {u_r}}{\partial {z}}u_\phi+\frac{\partial {u_\phi}}{\partial {z}}(r+u_r)\right) \\ \end{aligned}

To summarize we have, including both first and second order terms,

\begin{aligned}\begin{aligned}{d\mathbf{l}'}^2 - d\mathbf{x}^2&=2 dr dr \left( \frac{1}{{2}}\frac{\partial {u_m}}{\partial {r}} \frac{\partial {u_m}}{\partial {r}} +\frac{\partial {u_r}}{\partial {r}} \right) \\ &+2 r^2 d\phi d\phi \left( \frac{1}{{2 r^2}}\frac{\partial {u_m}}{\partial {\phi}} \frac{\partial {u_m}}{\partial {\phi}} + \frac{1}{{2 r^2}} \left( u_r^2 + u_\phi^2 \right) + \frac{u_r}{r}-\frac{1}{{r}}\frac{\partial {u_r}}{\partial {\phi}}\frac{u_\phi}{r}+\frac{1}{{r}}\frac{\partial {u_\phi}}{\partial {\phi}}\left(1+\frac{u_r}{r}\right)\right) \\ &+2 dz dz \left( \frac{1}{{2}}\frac{\partial {u_m}}{\partial {z}} \frac{\partial {u_m}}{\partial {z}} +\frac{\partial {u_z}}{\partial {z}} \right) \\ &+2 dr r d\phi \left( \frac{\partial {u_m}}{\partial {r}} \frac{1}{{r}}\frac{\partial {u_m}}{\partial {\phi}} +\frac{1}{{r}}\frac{\partial {u_r}}{\partial {\phi}} -\frac{u_\phi}{r}-\frac{\partial {u_r}}{\partial {r}}\frac{u_\phi}{r}+\frac{\partial {u_\phi}}{\partial {r}}\left(1 +\frac{u_r}{r}\right)\right) \\ &+2 r d\phi dz \left( \frac{\partial {u_m}}{\partial {z}} \frac{1}{{r}}\frac{\partial {u_m}}{\partial {\phi}} +\frac{1}{{r}}\frac{\partial {u_z}}{\partial {\phi}} -\frac{\partial {u_r}}{\partial {z}}\frac{u_\phi}{r}+\frac{\partial {u_\phi}}{\partial {z}}\left(1+\frac{u_r}{r}\right)\right) \\ &+2 dz dr \left( \frac{\partial {u_m}}{\partial {z}} \frac{\partial {u_m}}{\partial {r}} +\frac{\partial {u_r}}{\partial {z}} +\frac{\partial {u_z}}{\partial {r}} \right)\end{aligned}\end{aligned} \hspace{\stretch{1}}(3.60)

Factors of r have been pulled out so that the portions remaining in the braces are exactly the cylindrical tensor elements as given in the text (except also with the second order terms here). Observe that the pre-calculation of the general formula has allowed an on paper expansion of the cylindrical tensor without too much pain, and this time without requiring Mathematica.

Spherical tensor.



[1] L.D. Landau, EM Lifshitz, JB Sykes, WH Reid, and E.H. Dill. Theory of elasticity: Vol. 7 of course of theoretical physics. Physics Today, 13:44, 1960.

Posted in Math and Physics Learning. | Tagged: , , , , , , , | Leave a Comment »

A cylindrical Lienard-Wiechert potential calculation using multivector matrix products.

Posted by peeterjoot on May 1, 2011

[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.)]


A while ago I worked the problem of determining the equations of motion for a chain like object [1].
This was idealized as a set of N interconnected spherical pendulums. One of the aspects of that problem that I found fun was that it allowed me to use a new construct, factoring vectors into multivector matrix products, multiplied using the Geometric (Clifford) product. It seemed at the time that this made the problem tractable, whereas a traditional formulation was much less so. Later I realized that a very similar factorization was possible with matrices directly [2]. This was a bit disappointing since I was enamored by my new calculation tool, and realized that the problem could be tackled with much less learning cost if the same factorization technique was applied using plain old matrices.

I’ve now encountered a new use for this idea of factoring a vector into a product of multivector matrices. Namely, a calculation of the four vector Lienard-Wiechert potentials, given a general motion described in cylindrical coordinates. This I thought I’d try since we had a similar problem on our exam (with the motion of the charged particle additionally constrained to a circle).

The goal of the calculation.

Our problem is to calculate

\begin{aligned}A^0 &= \frac{q}{R^{*}} \\ \mathbf{A} &= \frac{q \mathbf{v}_c}{c R^{*}}\end{aligned} \hspace{\stretch{1}}(2.1)

where \mathbf{x}_c(t) is the location of the charged particle, \mathbf{r} is the point that the field is measured, and

\begin{aligned}R^{*} &= R - \frac{\mathbf{v}_c}{c} \cdot \mathbf{R} \\ R^2 &= \mathbf{R}^2 = c^2( t - t_r)^2 \\ \mathbf{R} &= \mathbf{r} - \mathbf{x}_c(t_r) \\ \mathbf{v}_c &= \frac{\partial {\mathbf{x}_c}}{\partial {t_r}}.\end{aligned} \hspace{\stretch{1}}(2.3)

Calculating the potentials for an arbitrary cylindrical motion.

Suppose that our charged particle has the trajectory

\begin{aligned}\mathbf{x}_c(t) = h(t) \mathbf{e}_3 + a(t) \mathbf{e}_1 e^{i \theta(t)}\end{aligned} \hspace{\stretch{1}}(3.7)

where i = \mathbf{e}_1 \mathbf{e}_2, and we measure the field at the point

\begin{aligned}\mathbf{r} = z \mathbf{e}_3 + \rho \mathbf{e}_1 e^{i \phi}\end{aligned} \hspace{\stretch{1}}(3.8)

The vector separation between the two is

\begin{aligned}\mathbf{R} &= \mathbf{r} - \mathbf{x}_c \\ &= (z - h) \mathbf{e}_3 + \mathbf{e}_1 ( \rho e^{i\phi} - a e^{i\theta} ) \\ &=\begin{bmatrix}\mathbf{e}_1 e^{i\phi} & - \mathbf{e}_1 e^{i\theta} & \mathbf{e}_3\end{bmatrix}\begin{bmatrix}\rho \\ a \\ z - h\end{bmatrix}\end{aligned}

Transposition does not change this at all, so the (squared) length of this vector difference is

\begin{aligned}\mathbf{R}^2 &=\begin{bmatrix}\rho &a & (z - h)\end{bmatrix}\begin{bmatrix}\mathbf{e}_1 e^{i\phi} \\ - \mathbf{e}_1 e^{i\theta} \\  \mathbf{e}_3\end{bmatrix}\begin{bmatrix}\mathbf{e}_1 e^{i\phi} & - \mathbf{e}_1 e^{i\theta} & \mathbf{e}_3\end{bmatrix}\begin{bmatrix}\rho \\ a \\ z - h\end{bmatrix} \\ &=\begin{bmatrix}\rho &a & (z - h)\end{bmatrix}\begin{bmatrix}\mathbf{e}_1 e^{i\phi} \mathbf{e}_1 e^{i\phi} & - \mathbf{e}_1 e^{i\phi} \mathbf{e}_1 e^{i\theta} & \mathbf{e}_1 e^{i\phi} \mathbf{e}_3 \\ - \mathbf{e}_1 e^{i\theta} \mathbf{e}_1 e^{i\phi} & \mathbf{e}_1 e^{i\theta} \mathbf{e}_1 e^{i\theta} & - \mathbf{e}_1 e^{i\theta} \mathbf{e}_3 \\  \mathbf{e}_3 \mathbf{e}_1 e^{i\phi} & -\mathbf{e}_3 \mathbf{e}_1 e^{i\theta} & \mathbf{e}_3 \mathbf{e}_3 \\ \end{bmatrix}\begin{bmatrix}\rho \\ a \\ z - h\end{bmatrix} \\ &=\begin{bmatrix}\rho &a & (z - h)\end{bmatrix}\begin{bmatrix}1 & - e^{i(\theta-\phi)} & \mathbf{e}_1 e^{i\phi} \mathbf{e}_3 \\ - e^{i(\phi -\theta)} & 1 & - \mathbf{e}_1 e^{i\theta} \mathbf{e}_3 \\  \mathbf{e}_3 \mathbf{e}_1 e^{i\phi} & -\mathbf{e}_3 \mathbf{e}_1 e^{i\theta} & 1 \\ \end{bmatrix}\begin{bmatrix}\rho \\ a \\ z - h\end{bmatrix} \\ \end{aligned}

A motivation for a Hermitian like transposition operation.

There are a few things of note about this matrix. One of which is that it is not symmetric. This is a consequence of the non-commutative nature of the vector products. What we do have is a Hermitian transpose like symmetry. Observe that terms like the (1,2) and the (2,1) elements of the matrix are equal after all the vector products are reversed.

Using tilde to denote this reversion, we have

\begin{aligned}(e^{i (\theta - \phi)})^{\tilde{}}&=\cos(\theta - \phi)+ (\mathbf{e}_1 \mathbf{e}_2)^{\tilde{}}\sin(\theta - \phi) \\ &=\cos(\theta - \phi)+ \mathbf{e}_2 \mathbf{e}_1\sin(\theta - \phi) \\ &=\cos(\theta - \phi)- \mathbf{e}_1 \mathbf{e}_2 \sin(\theta - \phi) \\ &=e^{-i (\theta -\phi)}.\end{aligned}

The fact that all the elements of this matrix, if non-scalar, have their reversed value in the transposed position, is sufficient to show that the end result is a scalar as expected. Consider a general quadratic form where the matrix has scalar and bivector grades as above, where there is reversion in all the transposed positions. That is

\begin{aligned}b^\text{T} A b\end{aligned} \hspace{\stretch{1}}(3.9)

where A = {\left\lVert{A_{ij}}\right\rVert}, a m \times m matrix where A_{ij} = \tilde{A_{ji}} and contains scalar and bivector grades, and b = {\left\lVert{b_i}\right\rVert}, a m\times 1 column matrix of scalars. Then the product is

\begin{aligned}\sum_{ij} b_i A_{ij} b_j&=\sum_{i<j} b_i A_{ij} b_j+\sum_{j<i} b_i A_{ij} b_j+\sum_{k} b_k A_{kk} b_k \\ &=\sum_{i<j} b_i A_{ij} b_j+\sum_{i<j} b_{j} A_{ji} b_i+\sum_{k} b_k A_{kk} b_k \\ &=\sum_{k} b_k A_{kk} b_k + 2 \sum_{i<j} b_i (A_{ij} + A_{ji}) b_j \\ &=\sum_{k} b_k A_{kk} b_k + 2 \sum_{i<j} b_i (A_{ij} + \tilde{A_{ij}}) b_j\end{aligned}

The quantity in braces A_{ij} + \tilde{A_{ij}} is a scalar since any of the bivector grades in A_{ij} cancel out. Consider a similar general product of a vector after the vector has been factored into a product of matrices of multivector elements

\begin{aligned}\mathbf{x} = \begin{bmatrix}a_1 & a_2 & \hdots & a_m\end{bmatrix}\begin{bmatrix}b_1 \\  b_2 \\  \dot{v}s \\  b_m\end{bmatrix}\end{aligned} \hspace{\stretch{1}}(3.10)

The (squared) length of the vector is

\begin{aligned}\mathbf{x}^2 &= (a_i b_i) (a_j b_j)  \\ &= (a_i b_i)^{\tilde{}} a_j b_j  \\ &= \tilde{b_i} \tilde{a_i} a_j b_j  \\ &= \tilde{b_i} (\tilde{a_i} a_j) b_j.\end{aligned}

It is clear that we want a transposition operation that includes reversal of its elements, so with a general factorization of a vector into matrices of multivectors \mathbf{x} = A b, it’s square will be \mathbf{x} = {\tilde{b}}^\text{T} {\tilde{A}}^\text{T} A b.

As with purely complex valued matrices, it is convenient to use the dagger notation, and define

\begin{aligned}A^\dagger = \tilde{A}^\text{T}\end{aligned} \hspace{\stretch{1}}(3.11)

where \tilde{A} contains the reversed elements of A. By extension, we can define dot and wedge products of vectors expressed as products of multivector matrices. Given \mathbf{x} = A b, a row vector and column vector product, and \mathbf{y} = C d, where each of the rows or columns has m elements, the dot and wedge products are

\begin{aligned}\mathbf{x} \cdot \mathbf{y} &= \left\langle{{ d^\dagger C^\dagger A b }}\right\rangle \\ \mathbf{x} \wedge \mathbf{y} &= {\left\langle{{ d^\dagger C^\dagger A b }}\right\rangle}_{2}.\end{aligned} \hspace{\stretch{1}}(3.12)

In particular, if b and d are matrices of scalars we have

\begin{aligned}\mathbf{x} \cdot \mathbf{y} &= d^\text{T} \left\langle{{C^\dagger A}}\right\rangle b = d^\text{T} \frac{C^\dagger A + A^\dagger C}{2} b \\ \mathbf{x} \wedge \mathbf{y} &= d^\text{T} {\left\langle{{C^\dagger A}}\right\rangle}_{2} b = d^\text{T} \frac{C^\dagger A - A^\dagger C}{2} b.\end{aligned} \hspace{\stretch{1}}(3.14)

The dot product is seen as a generator of symmetric matrices, and the wedge product a generator of purely antisymmetric matrices.

Back to the problem

Now, returning to the example above, where we want \mathbf{R}^2. We’ve seen that we can drop any bivector terms from the matrix, so that the squared length can be reduced as

\begin{aligned}\mathbf{R}^2 &=\begin{bmatrix}\rho &a & (z - h)\end{bmatrix}\begin{bmatrix}1 & - e^{i(\theta-\phi)} & 0 \\ - e^{i(\phi -\theta)} & 1 & 0 \\ 0 & 0 & 1 \\ \end{bmatrix}\begin{bmatrix}\rho \\ a \\ z - h\end{bmatrix} \\ &=\begin{bmatrix}\rho &a & (z - h)\end{bmatrix}\begin{bmatrix}1 & - \cos(\theta-\phi) & 0 \\ - \cos(\theta -\phi) & 1 & 0 \\ 0 & 0 & 1 \\ \end{bmatrix}\begin{bmatrix}\rho \\ a \\ z - h\end{bmatrix} \\ &=\begin{bmatrix}\rho &a & (z - h)\end{bmatrix}\begin{bmatrix}\rho - a \cos(\theta - \phi) \\ - \rho \cos(\theta - \phi) + a \\ z - h\end{bmatrix}\end{aligned}

So we have

\begin{aligned}\mathbf{R}^2 = \rho^2 + a^2 + (z -h)^2 - 2 a \rho \cos(\theta - \phi)\end{aligned} \hspace{\stretch{1}}(3.16)

Now consider the velocity of the charged particle. We can write this as

\begin{aligned}\frac{d \mathbf{x}_c}{dt} = \begin{bmatrix}\mathbf{e}_3 & \mathbf{e}_1 e^{i \theta} & \mathbf{e}_2 e^{i\theta}\end{bmatrix}\begin{bmatrix}\dot{h} \\ \dot{a} \\ a \dot{\theta}\end{bmatrix}\end{aligned} \hspace{\stretch{1}}(3.17)

To compute \mathbf{v}_c \cdot \mathbf{R} we have to extract scalar grades of the matrix product

\begin{aligned}\left\langle{{\begin{bmatrix}\mathbf{e}_1 e^{i\phi} \\ - \mathbf{e}_1 e^{i\theta} \\  \mathbf{e}_3\end{bmatrix}\begin{bmatrix}\mathbf{e}_3 & \mathbf{e}_1 e^{i \theta} & \mathbf{e}_2 e^{i\theta}\end{bmatrix}}}\right\rangle&=\left\langle{{\begin{bmatrix}\mathbf{e}_1 e^{i\phi} \\ - \mathbf{e}_1 e^{i\theta} \\  \mathbf{e}_3\end{bmatrix}\begin{bmatrix}\mathbf{e}_3 & \mathbf{e}_1 e^{i \theta} & \mathbf{e}_2 e^{i\theta}\end{bmatrix}}}\right\rangle \\ &=\left\langle{{\begin{bmatrix}\mathbf{e}_1 e^{i\phi} \mathbf{e}_3 & \mathbf{e}_1 e^{i\phi} \mathbf{e}_1 e^{i \theta}  & \mathbf{e}_1 e^{i\phi} \mathbf{e}_2 e^{i\theta} \\ - \mathbf{e}_1 e^{i\theta} \mathbf{e}_3 & - \mathbf{e}_1 e^{i\theta} \mathbf{e}_1 e^{i \theta}  & - \mathbf{e}_1 e^{i\theta} \mathbf{e}_2 e^{i\theta} \\  \mathbf{e}_3 \mathbf{e}_3 & \mathbf{e}_3 \mathbf{e}_1 e^{i \theta}  & \mathbf{e}_3 \mathbf{e}_2 e^{i\theta} \\ \end{bmatrix}}}\right\rangle \\ &= \begin{bmatrix}0 & \cos(\theta-\phi)  & - \sin(\theta - \phi) \\ 0 & - 1  & 0 \\ 1 & 0 & 0 \\ \end{bmatrix}.\end{aligned}

So the dot product is

\begin{aligned}\mathbf{R} \cdot \mathbf{v} &=\begin{bmatrix}\rho &a & (z - h)\end{bmatrix}\begin{bmatrix}0 & \cos(\theta-\phi)  & - \sin(\theta - \phi) \\ 0 & - 1  & 0 \\ 1 & 0 & 0 \\ \end{bmatrix}\begin{bmatrix}\dot{h} \\ \dot{a} \\ a \dot{\theta}\end{bmatrix} \\ &=\begin{bmatrix}\rho &a & (z - h)\end{bmatrix}\begin{bmatrix}\dot{a} \cos(\theta - \phi) - a \dot{\theta} \sin(\theta - \phi) \\ - \dot{a} \\ \dot{h}\end{bmatrix} \\ &=(z - h) \dot{h} - \dot{a} a + \rho \dot{a} \cos(\theta - \phi) - \rho a \dot{\theta} \sin(\theta - \phi) \end{aligned}

This is the last of what we needed for the potentials, so we have

\begin{aligned}A^0 &= \frac{q}{\sqrt{\rho^2 + a^2 + (z -h)^2 - 2 a \rho \cos(\theta - \phi)} -(z - h) \dot{h}/c + a \dot{a}/c + \rho \cos(\theta - \phi) \dot{a}/c - \rho a \sin(\theta - \phi) \dot{\theta}/c} \\ \mathbf{A} &= \frac{ \dot{h} \mathbf{e}_3 + (\dot{a} \mathbf{e}_1 + a \dot{\theta} \mathbf{e}_2) e^{i\theta} }{c} A^0,\end{aligned} \hspace{\stretch{1}}(3.18)

where all the time dependent terms in the potentials are evaluated at the retarded time t_r, defined implicitly by the messy relationship

\begin{aligned}c(t - t_r) = \sqrt{(\rho(t_r))^2 + (a(t_r))^2 + (z -h(t_r))^2 - 2 a(t_r) \rho \cos(\theta(t_r) - \phi)} .\end{aligned} \hspace{\stretch{1}}(3.20)

Doing this calculation with plain old cylindrical coordinates.

It’s worth trying this same calculation without any geometric algebra to contrast it. I’d expect that the same sort of factorization could also be performed. Let’s try it

\begin{aligned}\mathbf{x}_c &= \begin{bmatrix}a \cos\theta \\ a \sin\theta \\ h\end{bmatrix}\\ \mathbf{r} &= \begin{bmatrix}\rho \cos\phi \\ \rho \sin\phi \\ z\end{bmatrix}\end{aligned} \hspace{\stretch{1}}(4.21)

\begin{aligned}\mathbf{R} &= \mathbf{r} - \mathbf{x}_c \\ &= \begin{bmatrix}\rho \cos\phi - a \cos\theta \\ \rho \sin\phi - a \sin\theta \\ z - h\end{bmatrix} \\ &=\begin{bmatrix}\cos\phi & - \cos\theta & 0 \\ \sin\phi & - \sin\theta & 0 \\ 0 & 0 & 1 \end{bmatrix}\begin{bmatrix}\rho \\ a \\ z - h\end{bmatrix}\end{aligned}

So for \mathbf{R}^2 we really just need to multiply out two matrices

\begin{aligned}\begin{bmatrix}\cos\phi & \sin\phi & 0 \\ -\cos\theta & - \sin\theta & 0 \\ 0 & 0 & 1 \end{bmatrix}\begin{bmatrix}\cos\phi & - \cos\theta & 0 \\ \sin\phi & - \sin\theta & 0 \\ 0 & 0 & 1 \end{bmatrix}&=\begin{bmatrix}\cos^2\phi + \sin^2\phi & -(\cos\phi \cos\phi + \sin\phi \sin\theta) & 0 \\ -(\cos\phi \cos\theta + \sin\theta \sin\phi) & \cos^2\theta + \sin^2\theta & 0 \\ 0 & 0 & 1 \end{bmatrix} \\ &=\begin{bmatrix}1 & - \cos(\phi - \theta) & 0 \\ - \cos(\phi - \theta) & 1  & 0 \\ 0 & 0 & 1\end{bmatrix} \\ \end{aligned}

So for \mathbf{R}^2 we have

\begin{aligned}\mathbf{R}^2&=\begin{bmatrix}\rho & a & (z -h) \end{bmatrix}\begin{bmatrix}1 & - \cos(\phi - \theta) & 0 \\ - \cos(\phi - \theta) & 1  & 0 \\ 0 & 0 & 1\end{bmatrix} \begin{bmatrix}\rho \\ a \\ z - h\end{bmatrix} \\ &=\begin{bmatrix}\rho & a & (z -h) \end{bmatrix}\begin{bmatrix}\rho - a \cos(\phi - \theta) \\ -\rho \cos(\phi - \theta) + a \\ z - h\end{bmatrix} \\ &= (z - h)^2 + \rho^2 + a^2 - 2 a \rho \cos(\phi - \theta)\end{aligned}

We get the same result this way, as expected. The matrices of multivector products provide a small computational savings, since we don’t have to look up the \cos\phi \cos\phi + \sin\phi \sin\theta = \cos(\phi - \theta) identity, but other than that minor detail, we get the same result.

For the particle velocity we have

\begin{aligned}\mathbf{v}_c &= \begin{bmatrix}\dot{a} \cos\theta - a \dot{\theta} \sin\theta \\ \dot{a} \sin\theta + a \dot{\theta} \cos\theta \\ \dot{h} \end{bmatrix} \\ &=\begin{bmatrix}\cos\theta & - \sin\theta & 0 \\ \sin\theta & \cos\theta & 0 \\ 0 & 0 & 1\end{bmatrix} \begin{bmatrix}\dot{a} \\ a \dot{\theta}  \\ \dot{h} \end{bmatrix}\end{aligned}

So the dot product is

\begin{aligned}\mathbf{v}_c \cdot \mathbf{R} &=\begin{bmatrix}\dot{a} & a \dot{\theta}  & \dot{h} \end{bmatrix}\begin{bmatrix}\cos\theta & \sin\theta & 0 \\ -\sin\theta & \cos\theta & 0 \\ 0 & 0 & 1\end{bmatrix}\begin{bmatrix}\cos\phi & - \cos\theta & 0 \\ \sin\phi & - \sin\theta & 0 \\ 0 & 0 & 1 \end{bmatrix}\begin{bmatrix}\rho \\ a \\ z - h\end{bmatrix} \\ &=\begin{bmatrix}\dot{a} & a \dot{\theta}  & \dot{h} \end{bmatrix}\begin{bmatrix}\cos\theta \cos\phi + \sin\theta \sin\phi & -\cos^2 \theta - \sin^2 \theta & 0 \\ -\cos\phi \sin\theta + \cos\theta \sin\phi & 0 \\ 0 & 0 & 1\end{bmatrix}\begin{bmatrix}\rho \\ a \\ z - h\end{bmatrix} \\ &=\begin{bmatrix}\dot{a} & a \dot{\theta}  & \dot{h} \end{bmatrix}\begin{bmatrix}\cos(\phi - \theta) & -1 & 0 \\ \sin(\phi - \theta) & 0 & 0 \\ 0 & 0 & 1\end{bmatrix}\begin{bmatrix}\rho \\ a \\ z - h\end{bmatrix} \\ &=\dot{h}(z - h) - \dot{a} a + \rho \dot{a} \cos(\phi - \theta) + \rho a \dot{\theta} \sin(\phi - \theta)\end{aligned}

Reflecting on two the calculation methods.

With a learning curve to both Geometric Algebra, and overhead required for this new multivector matrix formalism, it is definitely not a clear winner as a calculation method. Having worked a couple examples now this way, the first being the N spherical pendulum problem, and now this potentials problem, I’ll keep my eye out for new opportunities. If nothing else this can be a useful private calculation tool, and the translation into more pedestrian matrix methods has been seen in both cases to not be too difficult.


[1] Peeter Joot. Spherical polar pendulum for one and multiple masses (Take II) [online].

[2] Peeter Joot. {Lagrangian and Euler-Lagrange equation evaluation for the spherical N-pendulum problem} [online].

Posted in Math and Physics Learning. | Tagged: , | Leave a Comment »

PHY450HS1: Relativistic electrodynamics: some exam reflection.

Posted by peeterjoot on April 28, 2011

[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.)]

Charged particle in a circle.

From the 2008 PHY353 exam, given a particle of charge q moving in a circle of radius a at constant angular frequency \omega.

\item Find the Lienard-Wiechert potentials for points on the z-axis.
\item Find the electric and magnetic fields at the center.

When I tried this I did it for points not just on the z-axis. It turns out that we also got this question on the exam (but stated slightly differently). Since I’ll not get to see my exam solution again, let’s work through this at a leisurely rate, and see if things look right. The problem as stated in this old practice exam is easier since it doesn’t say to calculate the fields from the four potentials, so there was nothing preventing one from just grinding away and plugging stuff into the Lienard-Wiechert equations for the fields (as I did when I tried it for practice).

The potentials.

Let’s set up our coordinate system in cylindrical coordinates. For the charged particle and the point that we measure the field, with i = \mathbf{e}_1 \mathbf{e}_2

\begin{aligned}\mathbf{x}(t) &= a \mathbf{e}_1 e^{i \omega t} \\ \mathbf{r} &= z \mathbf{e}_3 + \rho \mathbf{e}_1 e^{i \phi}\end{aligned} \hspace{\stretch{1}}(1.1)

Here I’m using the geometric product of vectors (if that’s unfamiliar then just substitute

\begin{aligned}\{\mathbf{e}_1, \mathbf{e}_2, \mathbf{e}_3\} \rightarrow \{\sigma_1, \sigma_2, \sigma_3\}\end{aligned} \hspace{\stretch{1}}(1.3)

We can do that since the Pauli matrices also have the same semantics (with a small difference since the geometric square of a unit vector is defined as the unit scalar, whereas the Pauli matrix square is the identity matrix). The semantics we require of this vector product are just \mathbf{e}_\alpha^2 = 1 and \mathbf{e}_\alpha \mathbf{e}_\beta = - \mathbf{e}_\beta \mathbf{e}_\alpha for any \alpha \ne \beta.

I’ll also be loose with notation and use \text{Real}(X) = \left\langle{{X}}\right\rangle to select the scalar part of a multivector (or with the Pauli matrices, the portion proportional to the identity matrix).

Our task is to compute the Lienard-Wiechert potentials. Those are

\begin{aligned}A^0 &= \frac{q}{R^{*}} \\ \mathbf{A} &= A^0 \frac{\mathbf{v}}{c},\end{aligned} \hspace{\stretch{1}}(1.4)


\begin{aligned}\mathbf{R} &= \mathbf{r} - \mathbf{x}(t_r) \\ R = {\left\lvert{\mathbf{R}}\right\rvert} &= c (t - t_r) \\ R^{*} &= R - \frac{\mathbf{v}}{c} \cdot \mathbf{R} \\ \mathbf{v} &= \frac{d\mathbf{x}}{dt_r}.\end{aligned} \hspace{\stretch{1}}(1.6)

We’ll need (eventually)

\begin{aligned}\mathbf{v} &= a \omega \mathbf{e}_2 e^{i \omega t_r} = a \omega ( -\sin \omega t_r, \cos\omega t_r, 0) \\ \dot{\mathbf{v}} &= -a \omega^2 \mathbf{e}_1 e^{i \omega t_r} = -a \omega^2 (\cos\omega t_r, \sin\omega t_r, 0)\end{aligned} \hspace{\stretch{1}}(1.10)

and also need our retarded distance vector

\begin{aligned}\mathbf{R} = z \mathbf{e}_3 + \mathbf{e}_1 (\rho e^{i \phi} - a e^{i \omega t_r} ),\end{aligned} \hspace{\stretch{1}}(1.12)

From this we have

\begin{aligned}R^2 &= z^2 + {\left\lvert{\mathbf{e}_1 (\rho e^{i \phi} - a e^{i \omega t_r} )}\right\rvert}^2 \\ &= z^2 + \rho^2 + a^2 - 2 \rho a (\mathbf{e}_1 \rho e^{i \phi}) \cdot (\mathbf{e}_1 e^{i \omega t_r}) \\ &= z^2 + \rho^2 + a^2 - 2 \rho a \text{Real}( e^{ i(\phi - \omega t_r) } ) \\ &= z^2 + \rho^2 + a^2 - 2 \rho a \cos(\phi - \omega t_r)\end{aligned}


\begin{aligned}R = \sqrt{z^2 + \rho^2 + a^2 - 2 \rho a \cos( \phi - \omega t_r ) }.\end{aligned} \hspace{\stretch{1}}(1.13)

Next we need

\begin{aligned}\mathbf{R} \cdot \mathbf{v}/c&= (z \mathbf{e}_3 + \mathbf{e}_1 (\rho e^{i \phi} - a e^{i \omega t_r} )) \cdot  \left(a \frac{\omega}{c} \mathbf{e}_2 e^{i \omega t_r} \right) \\ &=a \frac{\omega }{c}\text{Real}(i (\rho e^{-i \phi} - a e^{-i \omega t_r} ) e^{i \omega t_r} ) \\ &=a \frac{\omega }{c}\rho \text{Real}( i e^{-i \phi + i \omega t_r} ) \\ &=a \frac{\omega }{c}\rho \sin(\phi - \omega t_r)\end{aligned}

So we have

\begin{aligned}R^{*} = \sqrt{z^2 + \rho^2 + a^2 - 2 \rho a \cos( \phi - \omega t_r ) }-a \frac{\omega }{c} \rho \sin(\phi - \omega t_r)\end{aligned} \hspace{\stretch{1}}(1.14)

Writing k = \omega/c, and having a peek back at 1.4, our potentials are now solved for

\begin{aligned}\boxed{\begin{aligned}A^0 &= \frac{q}{\sqrt{z^2 + \rho^2 + a^2 - 2 \rho a \cos( \phi - k c t_r ) }} \\ \mathbf{A} &= A^0 a k ( -\sin k c t_r, \cos k c t_r, 0).\end{aligned}}\end{aligned} \hspace{\stretch{1}}(1.24)

The caveat is that t_r is only specified implicitly, according to

\begin{aligned}\boxed{c t_r = c t - \sqrt{z^2 + \rho^2 + a^2 - 2 \rho a \cos( \phi - k c t_r ) }.}\end{aligned} \hspace{\stretch{1}}(1.16)

There doesn’t appear to be much hope of solving for t_r explicitly in closed form.

General fields for this system.


\begin{aligned}\mathbf{R}^{*} = \mathbf{R} - \frac{\mathbf{v}}{c} R,\end{aligned} \hspace{\stretch{1}}(1.17)

the fields are

\begin{aligned}\boxed{\begin{aligned}\mathbf{E} &= q (1 - \mathbf{v}^2/c^2) \frac{\mathbf{R}^{*}}{{R^{*}}^3} + \frac{q}{{R^{*}}^3} \mathbf{R} \times (\mathbf{R}^{*} \times \dot{\mathbf{v}}/c^2) \\ \mathbf{B} &= \frac{\mathbf{R}}{R} \times \mathbf{E}.\end{aligned}}\end{aligned} \hspace{\stretch{1}}(1.18)

In there we have

\begin{aligned}1 - \mathbf{v}^2/c^2 = 1 - a^2 \frac{\omega^2}{c^2} = 1 - a^2 k^2\end{aligned} \hspace{\stretch{1}}(1.19)


\begin{aligned}\mathbf{R}^{*} &= z \mathbf{e}_3 + \mathbf{e}_1 (\rho e^{i \phi} - a e^{i k c t_r} )-a k \mathbf{e}_2 e^{i k c t_r} R \\ &= z \mathbf{e}_3 + \mathbf{e}_1 (\rho e^{i \phi} - a (1 - k R i) e^{i k c t_r} )\end{aligned}

Writing this out in coordinates isn’t particularly illuminating, but can be done for completeness without too much trouble

\begin{aligned}\mathbf{R}^{*} = ( \rho \cos\phi - a \cos t_r + a k R \sin t_r,  \rho \sin\phi - a \sin t_r - a k R \cos t_r,  z )\end{aligned} \hspace{\stretch{1}}(1.20)

In one sense the problem could be considered solved, since we have all the pieces of the puzzle. The outstanding question is whether or not the resulting mess can be simplified at all. Let’s see if the cross product reduces at all. Using

\begin{aligned}\mathbf{R} \times (\mathbf{R}^{*} \times \dot{\mathbf{v}}/c^2) =\mathbf{R}^{*} (\mathbf{R} \cdot \dot{\mathbf{v}}/c^2) - \frac{\dot{\mathbf{v}}}{c^2}(\mathbf{R} \cdot \mathbf{R}^{*})\end{aligned} \hspace{\stretch{1}}(1.21)

Perhaps one or more of these dot products can be simplified? One of them does reduce nicely

\begin{aligned}\mathbf{R}^{*} \cdot \mathbf{R} &= ( \mathbf{R} - R \mathbf{v}/c ) \cdot \mathbf{R}  \\ &= R^2 - (\mathbf{R} \cdot \mathbf{v}/c) R \\ &= R^2 - R a k \rho \sin(\phi - k c t_r) \\ &= R(R - a k \rho \sin(\phi - k c t_r))\end{aligned}

\begin{aligned}\mathbf{R} \cdot \dot{\mathbf{v}}/c^2&=\Bigl(z \mathbf{e}_3 + \mathbf{e}_1 (\rho e^{i \phi} - a e^{i \omega t_r} ) \Bigr) \cdot(-a k^2 \mathbf{e}_1 e^{i \omega t_r} )  \\ &=- a k^2 \left\langle{{\mathbf{e}_1 (\rho e^{i \phi} - a e^{i \omega t_r} ) \mathbf{e}_1 e^{i \omega t_r} )  }}\right\rangle \\ &=- a k^2 \left\langle{{(\rho e^{i \phi} - a e^{i \omega t_r} ) e^{-i \omega t_r} )  }}\right\rangle \\ &=- a k^2 \left\langle{{\rho e^{i \phi - i \omega t_r} - a }}\right\rangle \\ &=- a k^2 ( \rho \cos(\phi - k c t_r) - a )\end{aligned}

Putting this cross product back together we have

\begin{aligned}\mathbf{R} \times (\mathbf{R}^{*} \times \dot{\mathbf{v}}/c^2)&=a k^2 ( a -\rho \cos(\phi - k c t_r) ) \mathbf{R}^{*} +a k^2 \mathbf{e}_1 e^{i k c  t_r} R(R - a k \rho \sin(\phi - k c t_r)) \\ &=a k^2 ( a -\rho \cos(\phi - k c t_r) ) \Bigl(z \mathbf{e}_3 + \mathbf{e}_1 (\rho e^{i \phi} - a (1 - k R i) e^{i k c t_r} )\Bigr) \\ &\qquad +a k^2 R \mathbf{e}_1 e^{i k c  t_r} (R - a k \rho \sin(\phi - k c t_r)) \end{aligned}


\begin{aligned}\phi_r = \phi - k c t_r,\end{aligned} \hspace{\stretch{1}}(1.22)

this can be grouped into similar terms

\begin{aligned}\begin{aligned}\mathbf{R} \times (\mathbf{R}^{*} \times \dot{\mathbf{v}}/c^2)&=a k^2 (a - \rho \cos\phi_r) z \mathbf{e}_3 \\ &+ a k^2 \mathbf{e}_1(a - \rho \cos\phi_r) \rho e^{i\phi} \\ &+ a k^2 \mathbf{e}_1\left(-a (a - \rho \cos\phi_r) (1 - k R i)+ R(R - a k \rho \sin \phi_r)\right) e^{i k c t_r}\end{aligned}\end{aligned} \hspace{\stretch{1}}(1.23)

The electric field pieces can now be collected. Not expanding out the R^{*} from 1.14, this is

\begin{aligned}\begin{aligned}\mathbf{E} &= \frac{q}{(R^{*})^3} z \mathbf{e}_3\Bigl( 1 - a \rho k^2 \cos\phi_r \Bigr) \\ &+\frac{q}{(R^{*})^3} \rho\mathbf{e}_1 \Bigl(1 - a \rho k^2 \cos\phi_r \Bigr) e^{i\phi} \\ &+\frac{q}{(R^{*})^3} a \mathbf{e}_1\left(-\Bigl( 1 + a k^2 (a - \rho \cos\phi_r) \Bigr) (1 - k R i)(1 - a^2 k^2)+ k^2 R(R - a k \rho \sin \phi_r)\right) e^{i k c t_r}\end{aligned}\end{aligned} \hspace{\stretch{1}}(1.24)

Along the z-axis where \rho = 0 what do we have?

\begin{aligned}R = \sqrt{z^2 + a^2 } \end{aligned} \hspace{\stretch{1}}(1.25)

\begin{aligned}A^0 = \frac{q}{R} \end{aligned} \hspace{\stretch{1}}(1.26)

\begin{aligned}\mathbf{A} = A^0 a k \mathbf{e}_2 e^{i k c t_r } \end{aligned} \hspace{\stretch{1}}(1.27)

\begin{aligned}c t_r = c t - \sqrt{z^2 + a^2 } \end{aligned} \hspace{\stretch{1}}(1.28)

\begin{aligned}\begin{aligned}\mathbf{E} &= \frac{q}{R^3} z \mathbf{e}_3 \\ &+\frac{q}{R^3} a \mathbf{e}_1\left(-( 1 - a^4 k^4 ) (1 - k R i)+ k^2 R^2 \right) e^{i k c t_r} \end{aligned}\end{aligned} \hspace{\stretch{1}}(1.29)

\begin{aligned}\mathbf{B} = \frac{ z \mathbf{e}_3 - a \mathbf{e}_1 e^{i k c t_r}}{R} \times \mathbf{E}\end{aligned} \hspace{\stretch{1}}(1.30)

The magnetic term here looks like it can be reduced a bit.

An approximation near the center.

Unlike the old exam I did, where it didn’t specify that the potentials had to be used to calculate the fields, and the problem was reduced to one of algebraic manipulation, our exam explicitly asked for the potentials to be used to calculate the fields.

There was also the restriction to compute them near the center. Setting \rho = 0 so that we are looking only near the z-axis, we have

\begin{aligned}A^0 &= \frac{q}{\sqrt{z^2 + a^2}} \\ \mathbf{A} &= \frac{q a k \mathbf{e}_2 e^{i k c t_r} }{\sqrt{z^2 + a^2}} = \frac{q a k (-\sin k c t_r, \cos k c t_r, 0)}{\sqrt{z^2 + a^2}} \\ t_r &= t - R/c = t - \sqrt{z^2 + a^2}/c\end{aligned} \hspace{\stretch{1}}(1.31)

Now we are set to calculate the electric and magnetic fields directly from these. Observe that we have a spatial dependence in due to the t_r quantities and that will have an effect when we operate with the gradient.

In the exam I’d asked Simon (our TA) if this question was asking for the fields at the origin (ie: in the plane of the charge’s motion in the center) or along the z-axis. He said in the plane. That would simplify things, but perhaps too much since A^0 becomes constant (in my exam attempt I somehow fudged this to get what I wanted for the v = 0 case, but that must have been wrong, and was the result of rushed work).

Let’s now proceed with the field calculation from these potentials

\begin{aligned}\mathbf{E} &= - \boldsymbol{\nabla} A^0 - \frac{1}{{c}} \frac{\partial {\mathbf{A}}}{\partial {t}} \\ \mathbf{B} &= \boldsymbol{\nabla} \times \mathbf{A}.\end{aligned} \hspace{\stretch{1}}(1.34)

For the electric field we need

\begin{aligned}\boldsymbol{\nabla} A^0  &= q \mathbf{e}_3 \partial_z (z^2 + a^2)^{-1/2} \\ &= -q \mathbf{e}_3 \frac{z}{(\sqrt{z^2 + a^2})^3},\end{aligned}


\begin{aligned}\frac{1}{{c}} \frac{\partial {\mathbf{A}}}{\partial {t}} =\frac{q a k^2 \mathbf{e}_2 \mathbf{e}_1 \mathbf{e}_2 e^{i k c t_r} }{\sqrt{z^2 + a^2}}.\end{aligned} \hspace{\stretch{1}}(1.36)

Putting these together, our electric field near the z-axis is

\begin{aligned}\mathbf{E} = q \mathbf{e}_3 \frac{z}{(\sqrt{z^2 + a^2})^3}+\frac{q a k^2 \mathbf{e}_1 e^{i k c t_r} }{\sqrt{z^2 + a^2}}.\end{aligned} \hspace{\stretch{1}}(1.37)

(another mistake I made on the exam, since I somehow fooled myself into forcing what I knew had to be in the gradient term, despite having essentially a constant scalar potential (having taken z = 0)).

What do we get for the magnetic field. In that case we have

\begin{aligned}\boldsymbol{\nabla} \times \mathbf{A}(z)&=\mathbf{e}_\alpha \times \partial_\alpha \mathbf{A} \\ &=\mathbf{e}_3 \times \partial_z \frac{q a k \mathbf{e}_2 e^{i k c t_r} }{\sqrt{z^2 + a^2}}  \\ &=\mathbf{e}_3 \times (\mathbf{e}_2 e^{i  k c t_r} ) q a  k \frac{\partial {}}{\partial {z}} \frac{1}{{\sqrt{z^2 + a^2}}} +q a  k \frac{1}{{\sqrt{z^2 + a^2}}} \mathbf{e}_3 \times (\mathbf{e}_2 \partial_z e^{i  k c t_r} ) \\ &=-\mathbf{e}_3 \times (\mathbf{e}_2 e^{i  k c t_r} ) q a  k \frac{z}{(\sqrt{z^2 + a^2})^3} +q a  k \frac{1}{{\sqrt{z^2 + a^2}}} \mathbf{e}_3 \times \left( \mathbf{e}_2 \mathbf{e}_1 \mathbf{e}_2 k c e^{i  k c t_r} \partial_z ( t - \sqrt{z^a + a^2}/c ) \right) \\ &=-\mathbf{e}_3 \times (\mathbf{e}_2 e^{i  k c t_r} ) q a  k \frac{z}{(\sqrt{z^2 + a^2})^3} -q a  k^2 \frac{z}{z^2 + a^2} \mathbf{e}_3 \times \left( \mathbf{e}_1 k e^{i  k c t_r} \right) \\ &=-\frac{q a k z \mathbf{e}_3}{z^2 + a^2} \times \left( \frac{ \mathbf{e}_2 e^{i k c t_r} }{\sqrt{z^2 + a^2}} + k \mathbf{e}_1 e^{i k c t_r} \right)\end{aligned}

For the direction vectors in the cross products above we have

\begin{aligned}\mathbf{e}_3 \times (\mathbf{e}_2 e^{i \mu})&=\mathbf{e}_3 \times (\mathbf{e}_2 \cos\mu - \mathbf{e}_1 \sin\mu) \\ &=-\mathbf{e}_1 \cos\mu - \mathbf{e}_2 \sin\mu \\ &=-\mathbf{e}_1 e^{i \mu}\end{aligned}


\begin{aligned}\mathbf{e}_3 \times (\mathbf{e}_1 e^{i \mu})&=\mathbf{e}_3 \times (\mathbf{e}_1 \cos\mu + \mathbf{e}_2 \sin\mu) \\ &=\mathbf{e}_2 \cos\mu - \mathbf{e}_1 \sin\mu \\ &=\mathbf{e}_2 e^{i \mu}\end{aligned}

Putting everything, and summarizing results for the fields, we have

\begin{aligned}\mathbf{E} &= q \mathbf{e}_3 \frac{z}{(\sqrt{z^2 + a^2})^3}+\frac{q a k^2 \mathbf{e}_1 e^{i \omega t_r} }{\sqrt{z^2 + a^2}} \\ \mathbf{B} &= \frac{q a k z}{ z^2 + a^2} \left( \frac{\mathbf{e}_1}{\sqrt{z^2 + a^2}} - k \mathbf{e}_2 \right) e^{i \omega t_r}\end{aligned} \hspace{\stretch{1}}(1.38)

The electric field expression above compares well to 1.29. We have the Coulomb term and the radiation term. It is harder to compare the magnetic field to the exact result 1.30 since I did not expand that out.

FIXME: A question to consider. If all this worked should we not also get

\begin{aligned}\mathbf{B} \stackrel{?}{=}\frac{z \mathbf{e}_3 - \mathbf{e}_1 a e^{i \omega t_r}}{\sqrt{z^2 + a^2}} \times \mathbf{E}.\end{aligned} \hspace{\stretch{1}}(1.40)

However, if I do this check I get

\begin{aligned}\mathbf{B} =\frac{q a z}{z^2 + a^2} \left( \frac{1}{{z^2 + a^2}} + k^2 \right) \mathbf{e}_2 e^{i \omega t_r}.\end{aligned} \hspace{\stretch{1}}(1.41)

Collision of photon and electron.

I made a dumb error on the exam on this one. I setup the four momentum conservation statement, but then didn’t multiply out the cross terms properly. This led me to incorrectly assume that I had to try doing this the hard way (something akin to what I did on the midterm). Simon later told us in the tutorial the simple way, and that’s all we needed here too. Here’s the setup.

An electron at rest initially has four momentum

\begin{aligned}(m c, 0)\end{aligned} \hspace{\stretch{1}}(2.42)

where the incoming photon has four momentum

\begin{aligned}\left(\hbar \frac{\omega}{c}, \hbar \mathbf{k} \right)\end{aligned} \hspace{\stretch{1}}(2.43)

After the collision our electron has some velocity so its four momentum becomes (say)

\begin{aligned}\gamma (m c, m \mathbf{v}),\end{aligned} \hspace{\stretch{1}}(2.44)

and our new photon, going off on an angle \theta relative to \mathbf{k} has four momentum

\begin{aligned}\left(\hbar \frac{\omega'}{c}, \hbar \mathbf{k}' \right)\end{aligned} \hspace{\stretch{1}}(2.45)

Our conservation relationship is thus

\begin{aligned}(m c, 0) + \left(\hbar \frac{\omega}{c}, \hbar \mathbf{k} \right)=\gamma (m c, m \mathbf{v})+\left(\hbar \frac{\omega'}{c}, \hbar \mathbf{k}' \right)\end{aligned} \hspace{\stretch{1}}(2.46)

I squared both sides, but dropped my cross terms, which was just plain wrong, and costly for both time and effort on the exam. What I should have done was just

\begin{aligned}\gamma (m c, m \mathbf{v}) =(m c, 0) + \left(\hbar \frac{\omega}{c}, \hbar \mathbf{k} \right)-\left(\hbar \frac{\omega'}{c}, \hbar \mathbf{k}' \right),\end{aligned} \hspace{\stretch{1}}(2.47)

and then square this (really making contractions of the form p_i p^i). That gives (and this time keeping my cross terms)

\begin{aligned}(\gamma (m c, m \mathbf{v}) )^2 &= \gamma^2 m^2 (c^2 - \mathbf{v}^2) \\ &= m^2 c^2 \\ &=m^2 c^2 + 0 + 0+ 2 (m c, 0) \cdot \left(\hbar \frac{\omega}{c}, \hbar \mathbf{k} \right)- 2 (m c, 0) \cdot \left(\hbar \frac{\omega'}{c}, \hbar \mathbf{k}' \right)- 2 \cdot \left(\hbar \frac{\omega}{c}, \hbar \mathbf{k} \right)\cdot \left(\hbar \frac{\omega'}{c}, \hbar \mathbf{k}' \right) \\ &=m^2 c^2 + 2 m c \hbar \frac{\omega}{c} - 2 m c \hbar \frac{\omega'}{c}- 2\hbar^2 \left(\frac{\omega}{c} \frac{\omega'}{c}- \mathbf{k} \cdot \mathbf{k}'\right) \\ &=m^2 c^2 + 2 m c \hbar \frac{\omega}{c} - 2 m c \hbar \frac{\omega'}{c}- 2\hbar^2 \frac{\omega}{c} \frac{\omega'}{c} (1 - \cos\theta)\end{aligned}

Rearranging a bit we have

\begin{aligned}\omega' \left( m + \frac{\hbar \omega}{c^2} ( 1 - \cos\theta ) \right) = m \omega,\end{aligned} \hspace{\stretch{1}}(2.48)


\begin{aligned}\omega' = \frac{\omega}{1 + \frac{\hbar \omega}{m c^2} ( 1 - \cos\theta ) }\end{aligned} \hspace{\stretch{1}}(2.49)

Pion decay.

The problem above is very much like a midterm problem we had, so there was no justifiable excuse for messing up on it. That midterm problem was to consider the split of a pion at rest into a neutrino (massless) and a muon, and to calculate the energy of the muon. That one also follows the same pattern, a calculation of four momentum conservation, say

\begin{aligned}(m_\pi c, 0) = \hbar \frac{\omega}{c}(1, \hat{\mathbf{k}}) + ( \mathcal{E}_\mu/c, \mathbf{p}_\mu ).\end{aligned} \hspace{\stretch{1}}(3.50)

Here \omega is the frequency of the massless neutrino. The massless nature is encoded by a four momentum that squares to zero, which follows from (1, \hat{\mathbf{k}}) \cdot (1, \hat{\mathbf{k}}) = 1^2 - \hat{\mathbf{k}} \cdot \hat{\mathbf{k}} = 0.

When I did this problem on the midterm, I perversely put in a scattering angle, instead of recognizing that the particles must scatter at 180 degree directions since spatial momentum components must also be preserved. This and the combination of trying to work in spatial quantities led to a mess and I didn’t get the end result in anything that could be considered tidy.

The simple way to do this is to just rearrange to put the null vector on one side, and then square. This gives us

\begin{aligned}0 &=\left(\hbar \frac{\omega}{c}(1, \hat{\mathbf{k}}) \right) \cdot\left(\hbar \frac{\omega}{c}(1, \hat{\mathbf{k}}) \right) \\ &=\left( (m_\pi c, 0) - ( \mathcal{E}_\mu/c, \mathbf{p}_\mu ) \right) \cdot \left( (m_\pi c, 0) - ( \mathcal{E}_\mu/c, \mathbf{p}_\mu ) \right) \\ &={m_\pi}^2 c^2 + {m_\nu}^2 c^2 - 2 (m_\pi c, 0) \cdot ( \mathcal{E}_\mu/c, \mathbf{p}_\mu ) \\ &={m_\pi}^2 c^2 + {m_\nu}^2 c^2 - 2 m_\pi \mathcal{E}_\mu\end{aligned}

A final re-arrangement gives us the muon energy

\begin{aligned}\mathcal{E}_\mu = \frac{1}{{2}} \frac{ {m_\pi}^2 + {m_\nu}^2 }{m_\pi} c^2\end{aligned} \hspace{\stretch{1}}(3.51)

Posted in Math and Physics Learning. | Tagged: , , , , , , , , | Leave a Comment »