This month’s installment focuses on that aspect of plasma dynamics that is strongly analogous with mechanical vibration: plasma oscillations.  Plasma oscillations were first predicted by Paul Langmuir in the 1920s and later observed in vacuum tubes of the day.

These oscillations manifest themselves as cooperative movements in a plasma’s charge distribution that behave sinusoidalally.  The harmonic nature of these deviations from perfect neutrality is basically due to an electrostatic restoring force that, for small amplitudes around equilibrium, is well approximated as quadratic.  Thus the electrostatic restoring force gives rise to a harmonic potential.  The analog with mechanical oscillations is completed with the identification of the electron and ion centers-of-mass velocities with the macroscopic kinetic energy of the mass on the spring.  Plasma oscillations provide strong support for the basic assumptions of electromagnetic theory and a nice example of how mechanics and electrodynamics meet in a physically relevant situation.

The aim here is to derive the formulae for the plasma oscillations for two related cases and then to interpret the results in terms of mechanical vibration with the familiar spring mass system.

Imagine that there is a slab of plasma of length $$\ell_0$$, with cross-sectional area $$A$$ and an initial number density of ions and electrons of $$n_0$$.  Since the condition of quasi-neutrality is in effect, there is no net electric field in any part of the plasma.  Now imagine that due to an outside agency, all of the electrons are displaced by the very small amount of $$\Delta x$$ to the left, leaving a net positive charge on the right and net negative charge on the left, with the thickness of these two ‘non-neutral’ regions is $$\Delta x$$.

How does the plasma respond to this out of equilibrium condition? The first approximation one can make is that the electrons are free to move but the ions are fixed in space.  This corresponds to the physical approximation that the ions have an infinite mass – an approximation that is well supported by the ~2000 – 1 ratio between the mass of a proton and that of an electron.

The electrons will feel a strong restoring force (red arrow) and they will move accordingly.  Since the ions are fixed, they will feel an equal but opposite restoring force (blue arrow) but they will not be able to respond.

The electric field that the electrons and ions will feel, which is most easily obtained from Gauss’s law, is non-zero between the strips of positive and negative charge in an exactly analogous fashion as to the textbook problem on parallel plate capacitors.  The density of ions in the red region is $$n_0$$ and the volume is $$A \Delta x$$.   Since $$\Delta x$$ is small compared to the other dimensions, the electric field is closely approximated as being entirely in the $$x$$-direction and the total flux through a Gaussian surface whose boundaries coincide with the ion region is due strictly to the left and right faces.  The right face flux is exactly zero since the flux from the electrons exactly cancels the flux from the ions.  The left face flux, which is equal to the total flux, is then

\[ \Phi_E = E A \; , \]

which is proportional to the total charge enclosed

\[ \Phi_E = \frac{e n_0 A \Delta x}{\epsilon_0} \; .\]

Equating gives

\[ E= \frac{e n_0 }{\epsilon_0 } \Delta x \; ,\]

which is the usual result for the field between a parallel plate capacitor ($$n_0 \Delta x$$ is the surface charge density).

As the electrons were all displaced uniformly and are experiencing a uniform force, they will move together as a unit, that is to say cooperatively.  The mass of the unit is simply the mass density $$m_e n_0$$ times the volume consumed $$V = A \ell_0$$.  The center of mass of the slab is at $$\ell/2- \Delta x$$ with the corresponding acceleration $$\Delta {\ddot x}$$.  The force on the slab is

\[ F = q_{tot} E = -e n_0 V \frac{e n_0}{ \epsilon_0} \Delta x\]

and the equation of motion of the entire slab is given by:

\[ m_e n_0 V \frac{d^2}{dt^2} \Delta x = -e n_0 V \frac{e n_0} {\epsilon_0} \Delta x \; .\]

Dividing out the common factors gives

\[ \frac{d^2}{dt^2} \Delta x = -\frac{e^2 n_0}{m_e \epsilon_0} \Delta x \; .\]

From this equation we would conclude that, for small displacements, the electron slab would oscillate at the angular frequency

\[ \omega_{pe} = \sqrt{ \frac{e^2 n_0}{m_e \epsilon_0} } \; , \]

which is called the electron plasma frequency (note it is really an angular frequency).

Relaxing the assumption that the ions don’t move makes the analysis somewhat more complicated but only changes the result in a small but crucial way.  The location of both electron and ion centers-of-mass now matter and must be tracked separately.  The picture of their motion is now

Since both species are free to move, the thickness of the exposed electron or ion regions (red and blue) is not determined solely by the movement of their respective centers-of-mass (the ion region pictured above is much larger than $$\Delta x_i$$) .  The total thickness of the exposed regions depends on the relative motion between both centers-of-mass

\[ t = \Delta x_i – \Delta x_e \; .\]

An initial analysis may suggest that there are 8 separate cases: positive or negative variations for each species (2×2) as well as whether the ions are more displaced versus the electrons.  In fact, there are only two cases based on whether $$t$$ is positive (ion strip is on the right) or whether $$t$$ is negative (ion strip is to the left).  The following table shows the 8 variations for $$\ell = 1$$.

The parameters with the ‘min’ designate the leftmost extent of the electron slab ($$e_{min}$$) and the ion slab ($$i_{min}$$).  The parameter $$t_{low} = i_{min} – e_{min}$$.  A negative value, such as in case 1, means that the ion slab has shifted more to the left than the electrons and that the exposed ion region (red) is now on the left while the exposed electron region (blue) is now on the right.  A positive value of $$t_{low}$$ means the converse.  The corresponding parameters at the rightmost side are defined in an analogous fashion.  The color coding assists in visualizing which species is in which location.  Note that the entire table can be summarized by $$\Delta x_i – \Delta x_e$$, whose magnitude and sign match $$t_{low}$$ and $$t_{high}$$.  Also note that the displacements in the table are purely geometric and that the actual dynamic displacements would not shift the total center-of-mass since the only forces in the problem are internal.

The electric field in between the two uncovered charge regions can also be obtained by Gauss’s law, but it is instructive to use the expression for the electric field due to a sheet of charge that and superposition. The electric field due to the ions will have a magnitude of

\[ E_i = \frac{e n_0}{2 \epsilon_0} \left(\Delta x_i – \Delta x_e \right) \]

using the standard result from elementary E&M.

Likewise the electric field due to the electrons will also have the same magnitude $$E_e = E_i$$. Their directions will be the same in the region between them (the neutral plasma) and opposite outside, so that a non-zero field with magnitude $$2 E_i$$ will result only in the middle region. The direction of the field depends on which one is on the right versus the left but it doesn’t matter for the sake of the dynamical analysis.
The equation of motion for each slab is then given by

\[ m_i \frac{d^2}{dt^2} \Delta x_i = -\frac{e^2 n_0}{\epsilon_0} \left(\Delta x_i – \Delta x_e\right) \]

and

\[ m_e \frac{d^2}{dt^2} \Delta x_e = \frac{e^2 n_0}{\epsilon_0} \left(\Delta x_i – \Delta x_e\right) \; . \]

Combining these two equations gives

\[ \frac{d^2}{dt^2} \left( \Delta x_i – \Delta x_e \right) = -\frac{e^2 n_0}{\epsilon_0}\left(\frac{1}{m_e} + \frac{1}{m_i} \right) \left(\Delta x_i – \Delta x_e\right) \]

from which we conclude that the frequency is now

\[ \omega_p^2 = \frac{e^2 n_0}{\epsilon_0}\left(\frac{1}{m_e} + \frac{1}{m_i} \right) = \frac{e^2 n_0}{\epsilon_0 m_e} + \frac{e^2 n_0}{\epsilon_0 m_i} \equiv \omega_{pe}^2 + \omega_{pi}^2 \; .\]

The above analysis has a direct analog in a mechanical system of two unequal masses $$m_1$$ and $$M_2$$ connected by a spring of spring constant $$k$$ and equilibrium length $$\ell_0$$.

Defining the deviations of each mass from the equilibrium positions as

\[ \Delta x_1 = x_1 – x_1^{(0)} \]

and

\[ \Delta x_2 = x_2 – x_2^{(0)} \]

the distance between them takes the form

\[ x_2 – x_1 = x_2^{(0)} + \Delta x_2 – x_1^{(0)} – \Delta x_1 \equiv \ell \; .\]

Combining these two definitions gives

\[ x_2 – x_1 = \Delta x_2 – \Delta x_1 + \ell_0 \equiv \delta \ell + \ell_0 \]

from which the potential follows as

\[ V = \frac{1}{2} k (\ell – \ell_0)^2 = \frac{1}{2} k (\Delta x_2 – \Delta x_1)^2 \; .\]

The Lagrangian takes the form

\[ L = \frac{1}{2} m_1 \Delta \dot x_1^2 + \frac{1}{2} M_2 \Delta \dot x_2^2 – \frac{1}{2} k (\Delta x_2 – \Delta x_1)^2 \; , \]

with the equations of motion becoming

\[ \frac{d}{dt} \frac{\partial L}{\partial \dot x_1} – \frac{\partial L}{\partial x_1} = m_1 \Delta \ddot x_1 + k(\Delta x_2 – \Delta x_1 ) = 0 \]

and

\[ \frac{d}{dt} \frac{\partial L}{\partial \dot x_2} – \frac{\partial L}{\partial x_2} = M_2 \Delta \ddot x_2 – k(\Delta x_2 – \Delta x_1 ) = 0 \]

These equations combine nicely into one expression for the dynamics of the change in the separation between the two masses

\[ \delta \ddot \ell = -\left( \frac{k}{m_1} + \frac{k}{M_2} \right) \delta \ell \; , \]

which immediately tells us that the frequency of the oscillations is given by

\[ \omega^2 = k \left(\frac{1}{m_1} + \frac{1}{M_2} \right) = \omega_1^2 + \omega_2^2 \; .\]

In the limit as $$ M_2 \rightarrow \infty $$, the equation of motion simplifies to

\[ \delta \ddot \ell + \frac{k}{m_1} \delta \ell = 0 \; ,\]

with frequency of

\[ \omega^2 = \frac{k}{m_1} = \omega_1^2 \; .\]

Thus there is a perfect analogy between the two pictures and the cooperative motion of the electrons and ions within the plasma can be interpreted in strictly mechanical terms.

Starting with this installment, the focus shifts from the collective phenomenon of waves in cold plasma to the analysis of single particle motion of a charged particle under the influence of given maagnetic and electric fields. In this column, the motion in a fixed magnetic field will be analyzed. This is the most basic type of motion for plasmas and is almost always treated in the early chapters of plasma textbooks. Oddly enough, it is almost always presented in a cumbersome fashion and the cited ‘exact analytic’ solutions are usually incomplete or inaccurate. The major reason for this situation seems to be the fact that the analytic solutions are of little use in plasma physics. Nonetheless, from a pedagogical point-of-view this disconnect is a problem that should be addressed. This column is devoted to just such an aim.

The starting point is the Lorentz force law from which we obtain the equation of motion

\[ m \frac{d}{dt} \vec v = q \left( \vec E + \vec v \times \vec B \right) \; ,\]

subject to the initial conditions $$\vec r(t=0) \equiv \vec r_0$$ and $$\vec v(t=0) \equiv \vec v_0$$.

The textbook first case assumes that $$\vec E = 0$$, $$\vec B = B \hat z$$, and that the initial conditions are non-zero only in the $$x-y$$ plane. Plugging all these assumptions in yields the two equations of motion
\[ {\dot v}_x = \omega_g v_y \]

and

\[ {\dot v}_y = – \omega_g v_x \; ,\]

where the gyrofrequency is defined as

\[ \omega_g = \frac{q B}{m} \; .\]

Most textbook then try to solve these equations in rather clumsy fashion or they pass over the solution and cite analytic solutions that are incomplete or inaccurate. These will be discussed in detail below.

A much better way to solve these equations is to actually employ the standard mathematical way of solving coupled initial value problems using the fundamental matrix. The technique starts with a rewriting of the equations of motion in matrix form

\[ \frac{d}{dt} \left[ \begin{array}{c} v_x \\ v_y \end{array} \right] = \omega_g \left[ \begin{array}{cc} 0 & 1 \\ -1 & 0 \end{array} \right] \left[ \begin{array}{c} v_x \\ v_y \end{array} \right] \; .\]

The formal solution is given by
\[ \left[ \begin{array}{c} v_x \\ v_y \end{array} \right] (t) = e^{\omega_g A t} \left[ \begin{array}{c} v_{x0} \\ v_{y0} \end{array} \right] \; , \]

with the matrix

\[ A = \left[ \begin{array}{cc} 0 & 1 \\ -1 & 0 \end{array} \right] \; .\]

The matrix exponetial term is easily computed once realizes that $$A$$ squares to

\[ A^2 = \left[ \begin{array}{cc} -1 & 0 \\ 0 & -1 \end{array} \right] \; , \]

which is proportional to the unit matrix.

This is convenient since the power series form of the matrix exponential separates into two terms, each of which converges to a trigonometric function

\[ e^{\omega_g A t} = \cos (\omega_g t) 1 + \sin (\omega_g t) A \; ,\]

where $$1$$ is the $$2 \times 2$$ unit matrix.

Expansion of the fundamental matrix in terms of components results in

\[ e^{\omega_g A t} = \left[ \begin{array}{cc} \cos(\omega_g t) & \sin(\omega_g t) \\ -\sin(\omega_g t) & \cos(\omega_g t) \end{array} \right] \; .\]

Solutions to the initial value problem, in terms of initial conditions, are

\[ v_x(t) = v_{x0} \cos(\omega_g t) + v_{y0} \sin(\omega_g t) \]

and

\[ v_y(t) = -v_{x0} \sin(\omega_g t) + v_{y0} \cos(\omega_g t) \; .\]

To get the particle trajectories, one must then integrate these expressions with respect to time. There are two different ways of doing this and the choice between them is governed mostly by taste.

The first way is to perform a definite integration of both sides between $$t_0$$ (taken, for convenience, to be zero) and current time $$t$$:

\[ \int_{0}^{t} d x = \int_{0}^{t} \; dt’ \; \left( v_{x0} \cos(\omega_g t’) + v_{y0} \sin(\omega_g t’) \right) \; .\]

Explicitly computing the integrals gives

\[ x(t) – x_0 = \left. \frac{v_{x0}}{\omega_g} \sin(\omega_g t’) \right|_{0}^{t} – \left. \frac{v_{y0}}{\omega_g} \cos(\omega_g t’) \right|_{0}{t} \; ,\]

which simplifies to

\[ x(t) = x_0 + \frac{v_{y0}}{\omega_g} + \frac{v_{x0}}{\omega_g} \sin(\omega_g t) – \frac{v_{y0}}{\omega_g} \cos(\omega_g t) \; .\]

In the second method, the integrals are treated as indefinite and a single constant of integration results, giving

\[ x(t) – A = \frac{v_{x0}}{\omega_g} \sin(\omega_g t) – \frac{v_{y0}}{\omega_g} \cos(\omega_g t) \; .\]

The value of

\[ A = x_0 + \frac{v_{y0}}{\omega_g} \]

is then obtained from the initial conditions, resulting in exactly the same result.

For completeness, the result for the motion in the $$y$$ direction is

\[ y(t) = y_0 – \frac{v_{x0}}{\omega_g} + \frac{v_{x0}}{\omega_g} \cos(\omega_g t) + \frac{v_{y0}}{\omega_g} \sin(\omega_g t) \]

Interestingly, of the various textbooks examined in the field, none solve for the particle motion in the way above nor do they cite a functional form for the particle motion so derived. The most common way of tackling the equation of motion is exemplified by the treatment in Baumjohann and Treumann. They differentiate the equations once to decouple them thus arriving at

\[ {\ddot v}_x = -\omega_g^2 v_x \]

and

\[ {\ddot v}_y = -\omega_g^2 v_y \; .\]

This form has a two-fold disadvantage. First, the gyrofrequency $$\omega_g$$ appears squared here, thus losing all notion that particles of different signs in charge gyrate in different directions. Second, the second-order nature of these equations implies that there needs to be two constants of integration for the velocity evolution and a corresponding additional two for the particle motion, when, in fact, Newton’s laws require only four constants total. They fumble at with tis a bit in switching back and forth between the gyrofrequency having a sign and being the absolute value. In addition, they cite the solution to the particle equations as

\[ x – x_0 = r_g \sin( \omega_g t ) \]

and

\[ y – y_0 = r_g \cos( \omega_g t ) \; ,\]

with

\[ r_g \frac{v_{\perp}}{|\omega_g|} \; .\]

Chen tackles the equations similarly, deriving the velocity double-dot form before presenting the solution for the particle motion of

\[ x – x_0 = -i e^{-\omega_g t} \]

and

\[ y – y_0 = \pm e^{-\omega_g t} \; ,\]

which become, upon taking the real part,

\[ x – x_0 = r_g \sin( \omega_g t ) \]

and

\[ y – y_0 = \pm r_g \cos( \omega_g t ) \; .\]

He make no mention as to when to use the two different signs.

In order to see that the solution is derived in this column is correct, a numerical simulation of the equations of motion was implemented in Python, using the IPython/Jupyter framework with numpy and scipy extensions. The figures below show the initial conditions – initial position (green dot) and initial velocity (green line) – and the resulting motion from the numerical integration (blue line) with the analytic solution (red and black dot for positive and negative charge).

Positive Charge

Negative Charge

The agreement is excellent, showing the necessity of the ‘extra’ terms here derived. Also note that the sense of rotation is opposite for the negative charge compared to the positive one.

In tha last installment, the structure of plasma waves propagating along the local magnetic field direction were derived and analyzed. This column is mostly focued on the similar type of analysis for the case of waves propagating perpendicular to the local magnetic field. Along the way, a specific point that was glossed over the in the last analysis will be explored more fully. There are two reasons for this. The first is strictly the desire for these columns to be as complete as possible and a subsequent reading of the last piece as a refresher for this one showed a deficit. The second reason is that the analysis for the perpendicular case requires that point.

As before, the central relationship is the tangent form of the dispersion relation

\[ \tan^2 \theta = -\frac{P(n^2-R)(n^2 – L)}{(Sn^2 – RL)(n^2 – P)} \; , \]

where $$\theta$$ is the angle between the direction of wave propagation and the magnetic field, $$n$$ is the index of refraction, and the terms $$S$$, $$P$$, $$R$$, and $$L$$ are

\[ S = 1 – \sum_s \frac{\omega_{ps}^2}{\omega^2 – \omega_{cs}^2} \; ,\]

\[ P = 1 – \sum_s \frac{\omega_{ps}^2}{\omega^2} \; ,\]

\[ R = 1 – \sum_s \frac{\omega_{ps}^2}{\omega(\omega+\omega_{cs})} \; ,\]

and

\[ L = 1 – \sum_s \frac{\omega_{ps}^2}{\omega(\omega-\omega_{cs})} \; ,\]

respectively.

Perpendicular propagation is obtained by setting $$\theta = \pi/2$$. Since the tangent diverges at $$\pi/2$$, the appropriate dispersion relations result when the denominator goes to zero (i.e. the expression has a pole). The poles of the tangent form of the dispersion relation occur at

\[ n^2 = \frac{RL}{S} \]

and

\[ n^2 = P \; .\]

It is interesting to note that, unlike, the parallel propagation case, there are only two dispersion relations rather than three. Finding the eigenvectors leads to a subtle point that was glossed over in the last column. The intent here is to perform the computation in the proper amount of detail and then to explain why the gloss was permissible in the case of parallel propagation.

To determine the corresponding eigenvectors, it is conceptually most clean to go back to the defining equation for the dispersion tensor, which is given by:

\[ \left( \left( {\vec n}^{\times} \right)^2 + \overleftrightarrow{K} \right) \left[ \begin{array}{c} E_x \\ E_y \\ E_z \end{array} \right] = 0 \; , \]

with

\[ \overleftrightarrow{K} = \left[ \begin{array}{ccc} S & -i D & 0 \\ i D & S &0 \\ 0 & 0 & P \end{array} \right] \; ,\]

and

\[\left( {\vec n}^{\times} \right)^2 = \left[ \begin{array}{ccc} -n^2 \cos^2 \theta & 0 & n^2 \sin \theta \cos \theta \\0 & -n^2 & 0 \\n^2 \sin \theta \cos \theta & 0 & -n^2 \sin ^2 \theta \end{array} \right] \; .\]

When the direction of propagation is perpendicular, then $$\theta = \pi/2$$, and the first term becomes

\[\left( {\vec n}^{\times} \right)^2 = \left[ \begin{array}{ccc} 0 & 0 & 0 \\0 & -n^2 & 0 \\0 & 0 & -n^2 \end{array} \right] \; .\]

The resulting three equations for the electric field are:

\[S E_x – i D E_y = 0 \; ,\]

\[i D E_x + S E_y = n^2 E_y \; , \]

and

\[ P E_z = n^2 E_z \; . \]

The first equation can be solve to express $$E_x$$ in terms of $$E_y$$ as:

\[ E_x = \frac{i D}{S} E_y \, \]

independent of the actual dispersion relation. The solutions of the other two equations, however, depends intimately on the particulars of the dispersion relation.

For the case where $$n^2 = P$$, the $$y$$-equation becomes

\[ \left( S^2 – D^2 – PS \right) E_y = 0 \; , \]

for which the only solution is $$E_y = 0$$, since $$S^2 – D^2 – PS \neq 0$$. The third equation becomes

\[ P E_z = P E_z \; ,\]

which is satisfied with $$E_z$$ being assigned any value $$E_0$$. The corresponding eigenvector

\[ \left[ \begin{array}{c} 0 \\ 0 \\ P \end{array} \right] \]

is aligned along the magnetic field while its direction of propagation is in the plane perpendicular. This is known, according to Gurnett and Bhattacharjee as the ordinary mode. Note that its frequency can be arbitrary unlike the $$P$$ mode in the parallel case.

In contrast, in the case where $$n^2 = RL/S$$, the $$y$$-equation

\[ \left( S^2 – D^2 – RL \right) E_y = 0 \]

has a non-trivial solution since

\[ S^2 – D^2 = RL \]

is an identity. And so the equation is satisfied with $$E_y$$ being assigned any value $$E_0$$. The $$z$$-equation

\[ P E_z = RL/S E_z \]

can only be satisfied by $$E_z = 0$$. The corresponding eigenvector is

\[ \left[ \begin{array}{c} \frac{i D}{S} E_0 \\ E_0 \\ 0 \end{array} \right] \; . \]

This wave mode has its electric field in the $$x-y$$ plane, meaning that it has components neccessarily along and perpendicular to the direction of propagation. The parallel component yields an electrostatic piece while the perpendicular component yields an electromagnetic one. In addtion, its magnitude is no longer independent of the local magnetic field, as was the other cases examine up to this point. The dependence on magetic field strength comes from the $$S$$ and $$D$$ terms that scale $$E_x$$. This wave is termed extraordinary and leads to the hybrid resonances that are discussed at length in Section 4.4.2 of Gurnett and Bhattacharjee.

One last note before closing out this column. Last month, a seemingly simpler way was presented to find the eigenvectors for wave propagation parallel to the local magnetic field. The reason for the shortcut used was that when $$\theta = 0$$

\[\left( {\vec n}^{\times} \right)^2 = \left[ \begin{array}{ccc} -n^2 & 0 & 0 \\0 & -n^2 & 0 \\ 0 & 0 & 0 \end{array} \right] \; .\]

The basic equation

\[ \left( \left( {\vec n}^{\times} \right)^2 + \overleftrightarrow{K} \right) \left[ \begin{array}{c} E_x \\ E_y \\ E_z \end{array} \right] = 0 \; \]

becomes

\[ \left[ \begin{array}{ccc} S & -i D & 0 \\ i D & S &0 \\ 0 & 0 & P \end{array} \right] \left[ \begin{array}{c} E_x \\ E_y \\ E_z \end{array} \right] = \left[ \begin{array}{ccc} n^2 & 0 & 0 \\0 & n^2 & 0 \\ 0 & 0 & 0 \end{array} \right] \left[ \begin{array}{c} E_x \\ E_y \\ E_z \end{array} \right] \; ,\]

which is a basic eigenvalue problem for $$\overleftrightarrow{K}$$ and this is why the last column used a simple method for finding only its eigenvectors.

The last column established the basic equations describing waves in cold magnetized plasma and derived some of the general results that can be deduced without finding explicit solutions. This column begins the detailed examination of the explicit solutions that are supported.

The central relationship for this analysis is the tangent form of the dispersion relation

\[ \tan^2 \theta = -\frac{P(n^2-R)(n^2 – L)}{(Sn^2 – RL)(n^2 – P)} \; , \]

where $$\theta$$ is the angle between the direction of wave propagation and the magnetic field, $$n$$ is the index of refraction, and the terms $$S$$, $$P$$, $$R$$, and $$L$$ are

\[ S = 1 – \sum_s \frac{\omega_{ps}^2}{\omega^2 – \omega_{cs}^2} \; ,\]

\[ P = 1 – \sum_s \frac{\omega_{ps}^2}{\omega^2} \; ,\]

\[ R = 1 – \sum_s \frac{\omega_{ps}^2}{\omega(\omega+\omega_{cs})} \; ,\]

and

\[ L = 1 – \sum_s \frac{\omega_{ps}^2}{\omega(\omega-\omega_{cs})} \; ,\]

respectively.

Selecting a particular value of $$\theta$$ sets the propagation of the wave relative to the magnetic field and constrains the dispersion relation. There are three cases to examine: 1) parallel propagation, 2) perpendicular propagation, and 3) oblique propagation. This column will focus on parallel propagation.

In this case the wave moves along the magnetic field and $$\theta = 0$$. Since $$\tan \theta = 0$$, the propagation modes can be read off of the tangent form by simply finding the zeros of the right-hand side. Doing so yields

\[ P = 0 \; ,\]

\[ n^2 = R \; ,\]

and
\[ n^2 = L \; .\]

While the dispersion tensor

\[ \overleftrightarrow{D}(\vec n, \omega) = \left[ \begin{array}{ccc} S-n^2 & -i D & 0 \\ i D & S-n^2 &0 \\ 0 & 0 & P \end{array} \right] \; .\]

was useful for getting the characteristic equation, the matrix that gives the eigenvectors is the dielectric tensor (since $$n^2$$, in this case, is essentially the eigenvalue – at least for $$x$$ and $$y$$). The dielectric tensor’s form is

\[ \overleftrightarrow{K} = \left[ \begin{array}{ccc} S & -i D & 0 \\ i D & S &0 \\ 0 & 0 & P \end{array} \right] \; ,\]

where $$D = \frac{1}{2}( R – L )$$.

A simple analysis shows that the eigenvector corresponding to $$P=0$$ must have the form

\[ \vec E_{P} = \left[ \begin{array}{c} 0 \\ 0 \\ E_{P} \end{array} \right] \; ,\]

where $$E_P$$ is arbitrary.

This mode corresponds to a longitudinal wave at a frequency equal to $$\pm \omega_p$$, the combined plasma frequency of the system:

\[ P = 0 \Rightarrow \omega^2 = \sum_s \omega_{ps}^2 \; .\]

Since the wave is longitudinal it can’t be electromagnetic in origin and is, instead, associated with electrostatic fluctuations of the charge density along the magnetic field direction.

The eigenvectors for the other two modes require a bit more work but it is obvious that the corresponding eigenvectors can only have non-zero components along the $$x$$- and $$y$$-directions. Hence the structure of the wave is transverse, with the wave’s electric and magnetic fields perpendicular to the external magnetic field. Using wxMaxima, the eigenvector/eigenvalue pairs are given by:

\[ S + D = R: \left[ \begin{array}{c} 1 \\ i \\ 0 \end{array} \right] \]

and

\[ S – D = L: \left[ \begin{array}{c} 1 \\ -i \\ 0 \end{array} \right] \; .\]

The difference in the signs of the imaginary components correspond to the differences in handedness (hence the symbols $$R$$ and $$L$$ for right-handed and left-handed, respectively).

The usual convention for pulling out the handedness comes from multiplying the eigenvector by the time evolution $$e^{-i\omega t}$$ and taking the real part. For the $$R$$-mode, this process results in time dependence of

\[ E_x = E_0 \cos(\omega t) \]

and

\[ E_y = E_0 \sin(\omega t) \; ,\]

which shows that the electric field rotates in a counter-clockwise or right-handed fashion.

For the $$L$$-mode, the sign on the $$y$$-term causes the rotation to be clockwise or left-handed.

As Gurnett and Bhattacharjee point out, further proof that the modes are electromagnetic comes from noting that Fourier version of Faraday’s law gives a non-zero magnetic field

\[ \vec B = \vec k \times \vec E / \omega \; . \]

They further point out that since the wave is transverse then there are no charge fluctuations since Gauss’s law (in Fourier space) reads

\[ \rho = \epsilon_0 i \vec k \cdot \vec E = 0 \; . \]

Having specified the basic nature of these transverse waves, the final step is to say something about their specific dispersion relations.

There are three regimes to consider: 1) high frequency, 2) low frequency, and 3) intermediate frequencies.

High frequency

The high-frequency behavior of both modes is easy to determine. The limit of the corresponding dispersion relations as $$\omega \rightarrow \infty$$ is

\[ \lim_{\omega \rightarrow \infty} R = \lim_{\omega \rightarrow \infty} \left( 1 – \sum_s \frac{\omega_{ps}^2}{\omega(\omega+\omega_{cs})} \right) = 1 \; ,\]

and

\[ \lim_{\omega \rightarrow \infty} L = \lim_{\omega \rightarrow \infty} \left( 1 – \sum_s \frac{\omega_{ps}^2}{\omega(\omega-\omega_{cs})} \right) = 1 \; .\]

In this limit, both modes have $$n^2 = 1$$, which corresponds to free space propagation. Physically, this means that the plasma, composed of massive particles, is unable to respond at the rate the electromagnetic wave is oscillating at and the medium become transparent.

Low frequency

The low-frequency behavior of the two modes is a bit more subtle and are not tackled directly since the $$\omega(\omega \pm \omega_{cs})$$ term in the denominator is difficult to determine. Instead, the best approach is to take the limit of $$S$$ and $$D$$ and then to infer the limits for $$R$$ and $$L$$.

The limit of $$S$$ is

\[ \lim_{\omega \rightarrow 0} S = \lim_{\omega \rightarrow 0} \left( 1 – \sum_s \frac{\omega_{ps}^2}{\omega^2 – \omega_{cs}^2} \right) = 1 + \frac{\omega_{ps}^2}{\omega_{cs}^2} \; .\]

The limit of $$D$$

\[ \lim_{\omega \rightarrow 0} D = \lim_{\omega \rightarrow 0} \left( 1 – \sum_s \frac{\omega_{cs} \omega_{ps}^2}{\omega(\omega^2 – \omega_{cs}^2)} \right) \]

is a bit more involved. In this limit $$\omega^2 – \omega_{cs}^2 = -\omega_{cs}^2$$ and the term can be simplified to

\[ \lim_{\omega \rightarrow 0} \left( 1 + \frac{1}{\omega} \sum_s \frac{\omega_{ps}^2}{\omega_{cs}} \right) \; .\]

The next step is to examine the last term in detail. Expanding leads to

\[ \sum_s \frac{\omega_{ps}^2}{\omega_{cs}} = \sum_s \frac{n e^2}{m_s \epsilon_0} \frac{m_s}{q_s B} = \frac{1}{\epsilon_0 B} \sum_s n q_s = 0 \; , \]

where the requirement of neutrality is used in the last step.

Thus the low-frequency limit of $$D$$ is zero. From this result immediately follows that the low-frequency limits of $$R$$ and $$L$$ are the same as for $$S$$. This limit is fairly common and is defined as the Alfven index of refraction

\[ n_A^2 = 1 + \sum_s \frac{\omega_{ps}^2}{\omega_{cs}^2} \; .\]

Intermediate Frequencies

At intermediate frequencies the index of refraction departs from the low and high frequency limits and can do two very interesting things. First, near either the electron or cyclotron frequency, the $$R$$ or $$L$$ modes can enter into resonance with the electrons or the ions, respectively. Because the cyclotron frequency depends on the sign of the charge, the species interact with the two modes very differently. The right-hand mode can preferentially interact with the electrons near the electron cyclotron frequency. When this happens the index of refraction diverges to positive infinity as approached from below and diverges to negative infinity as approached from above. This last obseravation means that the wave cannot propagate until in the region from the negative regime to the cutoff frequency where $$n^2 > 0$$. Similar bevaior occurs for the $$L$$ mode with the ions with as many cutoffs and resonances as their are species.

The cutoff frequencies correspond to either $$R=0$$ or $$L=0$$. Since ion motion can be ignored at frequencies near the electron cycltron frequency, the $$R$$-mode cutoff is

\[ \omega_{R=0} = \frac{|\omega_{ce}|}{2} + \sqrt{\left(\frac{\omega_{ce}}{2}\right)^2 + \omega_{pe}^2} \; . \]

In the last installment, we examined a basic model of waves in cold plasmas. The two defining requirements were that the thermal fluctuations were small and ignorable (the cold piece) and that the plasma was not subjected to a magnetic field (unmagnetized). In this column, the requirement for a non-zero magnetic field ($$\vec B = 0$$) is relaxed while the requirement that the plasma is cold stays. The resulting types of wave motion are a lot richer. In contrast to the unmagnitized case, where there are two isotropic populations of electrons and ions, the presence of the magnetic field breaks the isotropy and motion of the charged particles is strongly different along the field lines as opposed to perpendicular to them.

Following again Gurnett and Bhattacharjee, the magnetic field takes on the form of

\[ \vec B = \vec B^{(0)} + \vec B^{(1)} \; ,\]

where $$\vec B^{(0)}$$ is a constant vector field and $$\vec B^{(1)}$$ carries all the spatial and temporal variations, assumed to be small compared to the constant piece. A similar form results for the density of each species. As in the unmagnetized case, the particle velocities and the electric field are zero to zeroth order.

The equations of motion for particles of the individual species, calculated from the Lorentz force law, is

\[ m_s \frac{\partial \vec v_s^{(1)}}{\partial t} = e_s \left[ \vec E^{(1)} + \vec v_s^{(1)} \times \vec B^{(0)} \right] \; .\]

Taking the Fourier transform of this equation leads to the three scalar equations
\[ – i \omega m_s v_{sx} = e_s(E_x + v_{sy}B ) \; ,\]
\[ – i \omega m_s v_{sy} = e_s(E_y – v_{sx}B ) \; ,\]
and
\[ – i \omega m_s v_{sz} = e_s E_z \; ,\]

in frequency space.

Dividing both sides by $$-i m_s$$ and using the cyclotron frequency

\[ \omega_{cs} = \frac{e_s B}{m_s} \]

on the right-hand side, allows for these three equations to be written in the matrix form as

\[ \left[ \begin{array}{ccc} -i \omega & – \omega_{cs} & 0 \\ \omega_{cs} & -i \omega & 0 \\ 0 & 0 & -i \omega \end{array} \right] \left[ \begin{array}{c} v_{sx} \\ v_{sy} \\ v_{sz} \end{array} \right] = \frac{e_s}{m_s} \left[ \begin{array}{c} E_x \\ E_y \\ E_z \end{array} \right] \; \]

or more compactly as

\[ {\mathbf M} \cdot \vec v = \vec E \; ,\]

where it is understood that all quantities are the Fourier Transforms of the corresponding state-space terms.

The goal is to express the velocties in terms of the electric field and, from them, the current density. To this end, the inverse of the matrix $${\mathbf M}$$ on the left-hand side is

\[{\mathbf M}^{-1} \equiv \left[ \begin{array}{ccc} \frac{i\omega}{\omega^2 – \omega_{cs}^2} & -\frac{\omega_{cs}}{\omega^2 – \omega_{cs}^2} & 0 \\ \frac{\omega_{cs}}{\omega^2 – \omega_{cs}^2} & \frac{i\omega}{\omega^2 – \omega_{cs}^2} & 0 \\0 & 0& \frac{i}{\omega} \end{array} \right] \; \]

from which the velocities can be written as

\[ \vec v = {\mathbf M}^{-1} \vec E \]

and the current densisty as

\[ \vec J = \sum_{s} n_s^{(0)}e_s \vec v \; .\]

Assuming the generalized Ohm’s law

\[ \vec J = \overleftrightarrow{\sigma} \vec E \; ,\]

immediately gives the form of the conductivity tensor

\[ \overleftrightarrow{\sigma} = \sum_s \frac{n_s^{(0)} e_s^2}{m_s} \left[ \begin{array}{ccc} \frac{i\omega}{\omega^2 – \omega_{cs}^2} & -\frac{\omega_{cs}}{\omega^2 – \omega_{cs}^2} & 0 \\ \frac{\omega_{cs}}{\omega^2 – \omega_{cs}^2} & \frac{i\omega}{\omega^2 – \omega_{cs}^2} & 0 \\0 & 0& \frac{i}{\omega} \end{array} \right] \; .\]

Plugging this last relation into the expression for the dielectric tensor

\[ \overleftrightarrow{K} = \overleftrightarrow{1} – \frac{\overleftrightarrow{\sigma}}{i \omega \epsilon_0} \]

gives an explicit form for the dielectric tensor for a magnetized, cold plasma

\[ \overleftrightarrow{K} = \left[ \begin{array}{ccc} 1 – \sum_s \frac{\omega_{ps}^2}{\omega^2 – \omega_{cs}^2} & -i \sum_s \frac{\omega_{ps}^2 \omega_{cs}}{\omega(\omega^2 – \omega_{cs}^2} & 0 \\ i \sum_s \frac{\omega_{ps}^2 \omega_{cs}}{\omega(\omega^2 – \omega_{cs}^2)}& 1 – \sum_s \frac{\omega_{ps}^2}{\omega^2 – \omega_{cs}^2} & 0\\0 &0 & 1 – \sum_s \frac{\omega_{ps}^2}{\omega^2} \end{array} \right] \; .\]

This expression, being unwieldy, suggests a shortcut notation be adopted. The usual convention is to define three quantities:

\[ S = 1 – \sum_s \frac{\omega_{ps}^2}{\omega^2 – \omega_{cs}^2} \; ,\]

\[ D = \sum_s \frac{\omega_{ps}^2 \omega_{cs}}{\omega(\omega^2 – \omega_{cs}^2)} \; ,\]

and

\[ P = 1 – \sum_s \frac{\omega_{ps}^2}{\omega^2} \; .\]

In terms of these expressions, the dielectric tensor is concisely written as

\[ \overleftrightarrow{K} = \left[ \begin{array}{ccc} S& -iD&0 \\iD & S&0 \\0 &0 &P \end{array} \right] \; .\]

To get the dispersion relation, we need to relate the wave propagation direction to the magnetic field direction. Generally, the propagation direction will not be parallel to the magnetic field but will be off by an angle $$\theta$$. The plane containing the two vectors can be taken, without loss of generality, to be the $$x-z$$ plane. In these terms, the propagation vector $$\vec n$$ takes the form

\[\vec n \equiv \left[ \begin{array}{c}n \sin \theta \\ 0 \\ n \cos \theta \end{array} \right] \; .\]

The ancilliary tensors associated with $$\vec n$$ take the form

\[ {\vec n}^{\times} = \left[ \begin{array}{ccc} 0 & n \cos \theta & 0\\-n \cos \theta & 0 & n \sin \theta \\ 0 &- n\sin \theta & 0 \end{array} \right] \]

and

\[ \left( {\vec n}^{\times} \right)^2 = \left[ \begin{array}{ccc} -n^2 \cos^2 \theta & 0 & n^2 \cos \theta \sin \theta \\ 0&-n^2 &0 \\ n^2 \cos \theta \sin \theta & 0 & -n^2 \sin^2 \theta \end{array} \right] \; .\]

Combining the last term with the dielectric tensor gives the form for the dispersion tensor

\[ \overleftrightarrow{D}(\vec n, \omega) = \left[ \begin{array}{ccc} S-n^2 \cos^2 \theta & -i D & n^2 \cos \theta \sin \theta \\ i D & S-n^2 &0 \\ n^2 \cos \theta \sin \theta & 0 & P – n^2 \sin^2 \theta \end{array} \right] \; .\]

We obtain the dispersion relations for the various waves by setting the determinant of $$\overleftrightarrow{D}(\vec n, \omega)$$ equal to zero. For reasons that will become obvious in the next treatment, the expressions $$S$$ and $$D$$ are related to the terms $$R$$ and $$L$$ (for right- and left-handed) through the transfomations

\[ R = S + D \]

and

\[ L = S – D \; .\]

It is a simple matter to invert these relations and express $$S$$ and $$D$$ in terms of $$R$$ and $$L$$ to give

\[ S = \frac{1}{2}[ R + L ] \]

and

\[ D = \frac{1}{2}[ R – L ] \; .\]

Of particular convenience is the further expression

\[ S^2 – D^2 = RL \; .\]

Now we are in position to calculate the determinant of the dispersion tensor. It is convenient to expand along the bottom row. Doing so yields

\[ |\overleftrightarrow{D}(\vec n, \omega)| = -n^4 \sin^2 \theta \cos^2 \theta (S-n^2) \\ + (P – n^2 \sin^2 \theta) \left[ (S-n^2)(S-n^2\cos^2 \theta) – D^2 \right]\; .\]

After some careful but straightforward expasions and simplifications and using $$S^2-D^2 = RL$$, the determinant simplifies to

\[ |\overleftrightarrow{D}(\vec n, \omega)| = [S\sin^2 \theta + P \cos^2 \theta ] n^4 \\ – [RL\sin^2 \theta +PS (1+\cos^2 \theta)] n^2 + RLP \; .\]

Note that terms proportional to $$n^6$$ have canceled out. As displayed, this expression is readily solved for $$n^2$$ using the quadratic formula with

\[ A = S\sin^2 \theta + P \cos^2 \theta \; , \]

\[ B = RL\sin^2 \theta +PS (1+\cos^2 \theta) \; ,\]

and

\[ C = RLP \; .\]

Some of the wave physics can be gleaned just by examining the discriminant

\[ F^2 = B^2 – 4AC \; . \]

The first term on the right-hand side can be expanded and simplified, particularly by using $$1+\cos^2 \theta = 2 – \sin^2 \theta$$, to

\[ B^2 = (RL – PS)^2 \sin^4 \theta + 4 P^2 S^2 \cos^2 \theta + 4RLPS \sin^2 \theta \; . \]

Likewise the last term becomes

\[ 4AC = 4RLPS \sin^2 \theta – 4RLP^2 \cos^2 \theta \; .\]

Combining, grouping, and simplifying yields

\[F^2 = B^2 – 4AC = (RL-PS)^2 \sin^4 \theta + 4P^2 D^2 \cos^2 \theta \; . \]

Clearly the discriminant is always positive-definite, which in turn implies that the roots are given by

\[ n^2 = \frac{-B \pm F}{2A} \; ,\]

which are always real quantities. Thus the index of refraction is either purely real or purely imaginary and either the waves propagate without damping or they don’t propagate at all.

The last result needed before mining for specific wave types comes from a simple but important manipulation of the characteristic equation by rewriting

\[ B = RL \sin^2 \theta + PS(\sin^2 \theta + 2 \cos^2 \theta) \]

and

\[ C = RLP(\sin^2 \theta + \cos^2 \theta) \; . \]

The steps for arriving at the desired result are these. Substitute the new expressions for $$B$$ and $$C$$ into the characteristic equation. Expand and collect on $$\sin^2 \theta$$ and $$\cos^2 \theta$$ to get

\[ \sin^2 \theta [Sn^4 – (RL+PS) n^2 +RLP] \\ + \cos^2\theta [Pn^4 -2 P S n^2 + RLP ] = 0 \; .\]

Solving for $$\tan^2 \theta$$ gives

\[ \tan^2 \theta = – \frac{[Pn^4 -2 P S n^2 + RLP ]}{[Sn^4 – (RL+PS) n^2 +RLP]} \; . \]

The denominator can be factored as

\[ Sn^4 – (RL + PS)n^2 + RLP = (Sn^2 – RL)(n^2 – P) \; .\]

The numerator is a bit trickier. To factor, first substitute $$2 S = R + L $$. This step leaves the numerator in the form

\[ Pn^4 – P(R+L)n^2 + RLP = P(n^4-(R+L)n^2 + RL) \; , \]

which leads to the factoring

\[ P(n^4 – (R+L)n^2 + RL ) = P(n^2 – R)(n^2 – L) \; .\]

Substituting these new expressions in yields the final, more convenient arrangement

\[ \tan^2 \theta = -\frac{P(n^2-R)(n^2 – L)}{(Sn^2 – RL)(n^2 – P)} \; . \]

In the next column, we start with this last relationship and derive various wave modes and their behavior.

Last column introduced the underlying formalism for describing waves in plasmas.  The key entities in the analysis are the dispersion tensor $$\overleftrightarrow{D}(\vec n, \omega)$$ and the dielectric tensor $$\overleftrightarrow{K}$$.  The latter quantity connects the behavior of the plasma to the underlying particle motions via the conductivity tensor $$\overleftrightarrow{\sigma}$$.  Setting the determinant of the dispersion tensor to zero results in the conditions for wave propagation; that is to say the dispersion relation that connects wavelength and frequency response.

The overall approach is to find a tractable, linear set of equations that related the state variables to each other.  Since one condition for a plasma is the maintenance of neutrality, at least at fairly large scales, the assumption will be that the state variables consist of a constant piece at zeroth order plus a first-order piece that acts as perturbations:  $$\vec B = \vec B^{(0)} + \vec B^{(1)} + \ldots$$.  The perturbation terms (e.g. $$\vec B^{(1)}$$) will then carry all of the action.

This installment looks at the most basic perturbation scenario available – waves in an unmagnetized cold plasma (see Section 4.3 of Gurnett and Bhattacharjee).  The technical conditions the plasma must satisfy to be classified this way are that its zeroth-order magnetic field $$\vec B^{(0)} = 0$$ (unmagnetized) and that its zeroth order particle velocities $$\vec v_s^{(0)}$$ and electric field $$\vec E^{(0)}$$ are also zero.

The particle density is also expanded perturbatively as

\[ n_s = n_s^{(0)} + n_s^{(1)} + \ldots \]

subject to the quasi-neutrality condition

\[ \sum_s e_s n_s^{(0)} = 0 \; .\]

The index $$s$$ specifies the species of the plasma: $$s=e$$ for electrons, $$s=i$$ for ions or $$s=h$$ for hydrogen, $$s=he$$ for helium, and so on.

As discussed last week, the approach is to start with the individual particle motions and derive a chain of results that link these with the conductivity tensor, then the dielectric tensor, then the dispersion tensor, and finally with the latter’s determinant.

For an unmagnetized cold plasma, the particle equations of motion (Lorenz force law) are:

\[ m_s \frac{d \vec v_s^{(1)}}{dt} = e_s \left[ \vec E^{(1)} + \vec v^{(1)} \times \vec B^{(1)} \right] \]

and

\[ \vec J = \sum_s e_s n_s \vec v_s = \sum_s e_s (n_s^{(0)} + n_s^{(1)}) \vec v_s^{(1)} \; . \]

Embracing the usual spirit of perturbative analysis, we throw away anything involving products of any first-order terms.  The resulting equations simplify considerably to become:

\[ m_s \frac{d \vec v_s^{(1)}}{dt} = e_s \vec E^{(1)} \]

and

\[ \vec J = \sum_s e_s n_s^{(0)} \vec v_s^{(1)} \; .\]

The Fourier transform of the differential equations yields:

\[ i \omega m_s \vec v_s^{(1)} = e_s \vec E^{(1)} \; .\]

Eliminating $$\vec v_s^{(1)}$$ from $$\vec J$$ gives

\[ \vec J = \sum_s e_s n_s^{(0)} i \frac{e_s}{\omega m_s} \vec E^{(1)} \; , \]

from which one immediately reads off the conductivity tensor as

\[ \overleftrightarrow{\sigma} = \sum_s \frac{e_s^2 n_s^{(0)} }{-i \omega m_s} \overleftrightarrow{1} \; . \]

The dielectric tensor follows immediately as

\[ \overleftrightarrow{K} = \overleftrightarrow{1}\left[1 – \sum_s \frac{e_s^2 n_s}{\omega^2 \epsilon_0 m_s} \right] \; .\]

Note that the perturbation order notation has been dropped for notational convenience.

The plasma frequency for an individual species $$\omega_{ps}$$ is a fundamental quantity that is defined as

\[ \omega_{ps}^2 = \frac{e_s^2 n_s}{\epsilon_0 m_s} \; . \]

The total plasma frequency of then defined as

\[ \omega_p^2 = \sum_s \omega_{ps}^2 \; .\]

Now, as there is no other physical direction in the problem, the direction of propagation can be taken, without loss of generality, to be along the $$z$$-axis. Thus the index of refraction vector is

\[ \vec n = \left[ 0, 0, n \right]^{T} \; .\]

The matrix $$(\vec n ^{\times})^2$$ simplifies to

\[ (\vec n ^{\times})^2 = \left[ \begin{array}{ccc} -n^2 & 0 & 0 \\ 0 & -n^2 & 0 \\ 0 & 0 & 0 \end{array} \right] \; . \]

The dispersion tensor is then

\[ \overleftrightarrow{D}(\vec n,\omega)= \left[ \begin{array}{ccc} 1 – \omega_p^2/\omega^2 -n^2 & 0 & 0 \\ 0 & 1 – \omega_p^2/\omega^2 -n^2 & 0 \\ 0 & 0 & 1 – \omega_p^2/\omega^2 \end{array} \right] \; . \]

Demanding that the determinant is zero gives the characteristic equation

\[ \left(1 – \frac{\omega_p^2}{\omega^2} – n^2 \right)^2 \left(1 – \frac{\omega_p^2}{\omega^2} \right) = 0 \; . \]

There are clearly two roots. The first, simple root, along the direction of propagation is called the longitudinal mode and has the dispersion relation

\[ \omega^2 = \omega_p^2 \; .\]

The second root, which is a double root, is associated with transverse mode and has the dispersion relation

\[ \omega^2 = \omega_p^2 + c^2 k^2 \; , \]

where $$c$$ is the speed of light and $$k$$ is the wave number.

Once the dispersion relation is obtained, the analysis of the resulting wave motion follows from looking at the phase and group velocities.

The dispersion relation for the longitudinal mode is independent of the wave number and the resulting motion, electrostatic in origin, is really an oscillation at the plasma frequency of the charge fluctuations. No wave actually propagates since the group velocity $$d \omega / d k$$ is zero.

The situation for the transverse mode is quite different. At frequencies large compared to the plasma frequency the wave behaves like a free-space wave with the dispersion relation $$\omega = \pm c k$$. As the frequency approaches the plasma frequency the nature of the propagation changes dramatically and propagation ceases entirely below the plasma frequency. The group velocity for the wave

\[ v_g = \frac{d \omega}{dk} = c \sqrt{1 – \frac{\omega_p^2}{\omega^2} } \]

clearly shows the cutoff as $$\omega \rightarrow \omega_p$$. It is interesting to note that while the group velocity is always lees than $$c$$, the phase velocity

\[ v_p = \frac{\omega}{k} = \frac{c}{\sqrt{1 – \frac{\omega_p^2}{\omega^2} }} \]

is always greater than $$c$$.

Next column will look at what happens when the assumption of $$\vec B = 0 $$ is relaxed.

This week is the start of a several columns devoted to studying wave propagation in plasmas. The wave behavior in a plasma are very rich due to the collective behavior of the plasma that results from its desire to keep charge neutrality. The presentation here is strongly influenced by and closely follows the approach presented in Introduction to Plasma Physics with Space and Laboratory Applications by Gurnett and Bhattacharjee.

Waves in a plasma represent the propagation over long distances of a complicated patterns of interaction between the charge particles that make up the plasma and the electromagnetic fields present. Thus the starting point is the Lorentz force law

\[ m_s \frac{d \vec v_s}{d t} = q_s (\vec E + \vec v_s \times \vec B ) \; , \]

where $$s$$ can be thought of as a multi-index labeling both particle and species, and
Maxwell’s equations.

The particular form of Maxwell’s equations – microscopic versus macroscopic – that is best to use is a bit of a thorny question. On one hand, the electrons and ions that comprise a plasma are free to move about over very large distances. On the other, they act together to maintain charge neutrality. Gurnett and Bhattacharjee seem to split the difference by analyzing both the microscopic set

\[ \begin{array}{ccc} \nabla \cdot \vec E & = & \rho/\epsilon_0 = (\rho_{free}+\rho_{pol})/\epsilon_0 \\ \nabla \cdot \vec B & = & 0 \\ \nabla \times \vec E & = & – \partial_t \vec B \\ \nabla \times \vec B & = & \mu_0 \vec J + \mu_0 \epsilon_0 \partial_t \vec E \end{array} \]

and the macroscopic set

\[ \begin{array}{ccc} \nabla \cdot \vec D & = & \rho_{free} \\ \nabla \cdot \vec B & = & 0 \\ \nabla \times \vec E & = & – \partial_t \vec B \\ \nabla \times \vec H & = & \vec J_{free} + \partial_t \vec D \end{array} \; \]

and then enforcing a sort of consistency relation between them. The central notion is that all of the charges are to be counted as bound charges so that the usual linear relation between the displacement vector, the electric field, and the polarization – $$\vec D = \epsilon_0 \vec E + \vec P$$ holds. Gurnett and Bhattacharjee are not particularly clear on how to treat the magnetic constitutive relation other than to say

This half-and-half approach seems to be related to the points raised in the discussion in Section 3.2 of Chen’s book, where he points out that the dependence of the magnetic moment of the individual particles is given by

\[ \mu = \frac{m v_{\perp}^2}{2 B} \; , \]

which makes the magnetization, $$\vec M$$, nonlinear in $$|\vec B|$$, thus complicating the purely-macroscopic approach.

In any event, the Fourier transform is the usual tool used to tackle wave phenomena, since it provides a transform framework that handles both spatial and time variations in a convenient way. The Fourier transform replaces all time derivatives with multiplications by frequency and replaces all spatial derivatives with multiplication by a wave number:

\[ \partial_t \rightarrow – i \omega \]

\[ \nabla \cdot \rightarrow \vec k \cdot \]

\[ \nabla \times \rightarrow \vec k \times \]

Using this translations, the microscopic Maxwell equations become

\[ \begin{array}{ccc} i \hat k \cdot \vec E & = & \rho/\epsilon_0 \\ i \hat k \cdot \vec B & = & 0 \\ i \hat k \times \vec E & = & i \omega \vec B \\ i \hat k \times \vec B & = & \mu_0 \vec J – i \omega \mu_0 \epsilon_0 \vec E \end{array} \; ,\]

while the macroscopic equations become

\[ \begin{array}{ccc} i \hat k \cdot \vec E & = & 0 \\ i \hat k \cdot \vec B & = & 0 \\ i \hat k \times \vec E & = & i \omega \vec B \\ i \hat k \times \vec B & = & – i \omega \mu_0 \vec D \end{array} \; .\]

It is understood that the vector quantities in both set of equations are the Fourier transforms and, thus, depend on frequency and wave number and not time and position. Also note that the assumptions that a plasma has only bound charges and no magnetization results in no free density or current terms in the macroscopic set.

To close the set of equations, two constitutive relations need to be adopted. The first is the generalized Ohm’s law that relates the current to the electric field

\[ \vec J = \overleftrightarrow{\sigma} \cdot \vec E \; \]

through the conductivity tensor $$\overleftrightarrow{\sigma}$$.

The second one is a generalized dielectric relation between the displacement vector and the electric field

\[ \vec D = \epsilon_0 \overleftrightarrow{K} \cdot \vec E \]

through the dielectric tensor $$\overleftrightarrow{K}$$.

A formal connection between the conductivity and dielectric tensors can be made by equating the right-hand sides of Ampere’s law in microscopic and macroscopic form to yield

\[ \mu_0 \overleftrightarrow{\sigma} \cdot \vec E – i \omega \epsilon_0 \mu_0 \vec E = – i \omega \mu_0 \epsilon_0 \overleftrightarrow{K} \cdot \vec E \; .\]

A simple re-arrangement of terms results in the dielectric tensor express as

\[ \overleftrightarrow{K} = \overleftrightarrow{1} – \frac{ \overleftrightarrow{\sigma} }{i \omega \epsilon_0 } \; .\]

With the connection of the conductivity and dielectric tensors made, the microscopic equation now take a back seat to the macroscopic ones. The wave equation results when substituting Ampere’s equation into Faraday’s. Doing so yields

\[ \vec k (\vec k \times \vec E) + \frac{\omega^2}{c^2} \overleftrightarrow{K} \cdot \vec E = 0 \; .\]

Defining an index of refraction vector as

\[ \vec n = \frac{c}{\omega} \vec k\]

and using the fundamental relation

\[ \mu_0 \epsilon_0 = \frac{1}{c^2} \]

gives the form

\[ \vec n \times (\vec n \times \vec E) + \overleftrightarrow{K} \cdot \vec E = 0 \; .\]

The only way for this equation to be satisfied with a non-trivial electric field is for the operator that acts on it to be singular. If the equation can be written as

\[ \overleftrightarrow{D}(\vec n, \omega) \cdot \vec E = \vec 0 \]

then the determinant of $$\overleftrightarrow{D}(\vec n, \omega)$$ must be zero.

Unfortunately, $$\overleftrightarrow{D}(\vec n, \omega)$$ hides its structure with the complicated division into two pieces, one of which involves a double cross-product. Fortunately, there are convenient tools from mechanics that can tame the first term and join that which is split.

The first step is to recognize that a vector cross product can be represented in terms of matrix multiplication

\[ \vec n \times \doteq \left[ \begin{array}{ccc} 0 & -n_z & n_y \\n_z & 0 & -n_x \\ -n_y & n_x & 0 \end{array} \right] \; ,\]

which is often denoted as

\[ \vec n \times \equiv \vec n^{\times} \; .\]

This step is very convenient since it puts the cross product into matrix form and it rids the need for the parentheses. The double product immediately results from a second application of the matrix $$\vec n^{\times}$$ giving

\[ \vec n \times (\hat n \times \vec E) \doteq \left[ \begin{array}{ccc} -n_y^2 – n_z^2 & n_x n_y & n_x n_z \\n_x n_y & -n_x^2 – n_z^2 & n_y n_z \\n_x n_z & n_y n_z & -n_x^2 -n_y^2 \end{array} \right] \; .\]

Of course, this same result falls out from the expansion of $$\vec n \times (\vec n \times \vec E) = \vec n (\vec n \cdot \vec E) – \vec E (\vec n \cdot \vec n) $$ using BAC-CAB.

The operator $$\overleftrightarrow{D}(\vec n, \omega)$$, which I will call the dispersion tensor, is then recognized as the sum

\[ \overleftrightarrow{D}(\vec n, \omega) = (\vec n^{\times})^2 + \overleftrightarrow{K} = (\vec n^{\times})^2 + \overleftrightarrow{1} – \frac{\overleftrightarrow{\sigma}}{i \omega \epsilon_0} \; .\]

Much of the relevant behavior of the waves come from mining relationships that come from the requirement that the determinant of the dispersion tensor is zero. This is the coupled-equation form of the requirement that the wave operator $$\nabla^2 – c^{-2} \partial_t^2$$ is zero when operating on the a valid waveform (e.g. $$f(x-ct) + g(x+ct)$$).

The strategy for finding the dispersion tensor is as follows. First compute the conductivity tensor using some model of the plasma. The possible ways of tackling this seem to be single particle motion from the Lorentz law, magnetohydrodynamic methods, and kinetic theory. Next get the dielectric tensor by plugging the conductivity tensor into the relation above. Finally, add this piece to $$({\vec n^{\times}})^2$$.

Next week I’ll examine the implementation of this strategy for a cold plasma.

This week’s column is relatively short but sweet. It covers the general operations that can be carried out on the Gaussian family of integrals. Since these integrals are relatively easy to compute and they find applications in numerous fields they are a cornerstone of many physical disciplines and appear in many models of basic systems.

The prototype integral upon which the whole family is built is the traditional

\[ I = \int_{-\infty}^{\infty} dx \; e^{-a x^2} \; .\]

Although the derivation of this integral is well known, I’ll cover the computation for the sake of completeness.

Solving for the value of $$I$$ is much easier when one squares both sides to get

\[ I^2 = \int_{-\infty}^{\infty} dx \, e^{-a x^2} \int_{-\infty}^{\infty} dy \, e^{-a y^2} \; . \]

The next step is to put the repeated integral together as a double integral and then to change coordinates from Cartesian to plane polar. Following this program, one soon arrives at

\[ I^2 = \int_{0}^{\infty} \int_{0}^{2\pi}\, dr \, d\theta \, r e^{-ar^2} \; .\]

The $$\theta$$ integral is easily done leaving

\[ I^2 = 2 \pi \int_{0}^{\infty} \, dr \, r e^{-a r^2} \; .\]

Letting $$s = a r^2$$ gives $$ds = 2\, a\, r\, dr$$ with limits $$s(r=0) = 0$$ and $$s(r=\infty) = \infty$$. Plugging these in yields

\[ I^2 = \frac{\pi}{a} \int_{0}^{\infty} ds \, e^{-s} \; ,\]

which is a standard integral that can be immediate solved to give

\[ I^2 = \frac{\pi}{a} \left. e^{-s} \right|_{0}^{\infty} = \frac{\pi}{a} \]

or, more simply

\[ I = \int_{-\infty}^{\infty} \, dx \, e^{-a x^2} = \sqrt{ \frac{\pi}{a} } \; .\]

This can be generalized without much more work to the more fearsome-looking (but actually not much harder) integral family

\[ I_n = \int_{-\infty}^{\infty} \, dx \, x^n e^{-a x^2 + b x } \; n = 0, 1, 2, \ldots \; .\]

The first step is to evaluate $$I_0$$, which is the only true integral that will have to be done. The integral

\[ I_0 = \int_{-\infty}^{\infty} \, dx \, e^{-a x^2 + b x} \]

is tackled by completing the square:

\[ \int_{-\infty}^{\infty}\, dx e^{-a x^2 + b x} = \int_{-\infty}^{\infty} dx \; e^{-(\sqrt{a} x – \frac{b}{2\sqrt{a}})^2} e^{b^2/4a} \; . \]

The last term is a constant that can be brought out of the integral. The remaining piece is best computed by first defining

\[ y = \sqrt{a} x – \frac{b}{2\sqrt{a}} \; \]

and

\[ dy = \sqrt{a} dx \; .\]

The limits of integration stay the same the new form of the integral is

\[ I_0 = e^{b^2/4a} \frac{1}{\sqrt{a}} \int_{-\infty}^{\infty} \, dy e^{-y^2} = \sqrt{\frac{\pi}{a}} e^{b^2/4a} \; . \]

This result is so important it deserves its own box

All the other integrals come from the clever trick that says that powers of $$x$$ can be obtained from $$I_0$$ by differentiating with respect to either $$b$$ (all powers) or $$a$$ (even powers only).

For example, by definition

\[ I_1 = \int_{-\infty}^{\infty} \, dx \, x e^{-a x^2 + bx} \; ,\]

but, by construction

\[ \frac{\partial}{\partial b} I_0 = \frac{\partial}{\partial b} \int_{-\infty}^{\infty} \, dx \, e^{-a x^2 + bx} = \int_{-\infty}^{\infty} \, dx \, \frac{\partial}{\partial b} e^{-a x^2 + bx} \\ = \int_{-\infty}^{\infty} \, dx \, x e^{-a x^2 + bx} \; .\]

So

\[ I_1 = \frac{\partial}{\partial b} I_0 = \frac{\partial}{\partial b} \sqrt{\frac{\pi}{a}} e^{b^2/4a} \]

and the computation of this integral is done by performing a differentiation – a much easier task – to get

\[ I_1 = \frac{b}{4a} \sqrt{\frac{\pi}{a}} e^{b^2/4a} \; . \]

This generalizes immediately to

\[ I_n = \frac{\partial^n}{\partial b^n} I_0 \; .\]

Sometimes the integrals we wish to compute have $$b=0$$. In this case, which is fairly common, we can get the desired results almost as easily.

For $$n$$ odd, the integrand will be a product of an even function $$\exp({-ax^2})$$ and an odd function $$x^n$$ over an even domain. The result is exactly zero.

The $$n$$ even case comes from a similar differentiation with respect to $$a$$. For example,

\[ I_2 = -\frac{\partial}{\partial a} I_0 = \int_{-\infty}^{\infty} \, dx \, x^2 e^{-a x^2} \; .\]

Evaluating the derivative gives

\[ I_2 = -\frac{\partial}{\partial a} \sqrt{\frac{\pi}{a}} = \frac{1}{2} \sqrt{\frac{\pi}{a^3}} \; .\]

There are two consistency checks that are worth noting.

First, the recipe involving $$b$$ can still be used in the $$b = 0$$ case by taking the limit as $$b$$ approaches zero as the last step. In this limit,

\[ \lim_{b\rightarrow 0} I_1 = \lim_{b\rightarrow 0} \frac{b}{4a} \sqrt{\pi}{a} e^{b^2/4a} = 0 \; , \]

which checks with the partity argument given above.

Second, one would like to see that

\[ -\frac{\partial}{\partial a} I_0 = \frac{\partial^2}{\partial b^2} I_0 \]

for the computation of $$I_2$$. The left-hand side gives

\[ -\frac{\partial}{\partial a} I_0 = \frac{1}{2} \sqrt{\frac{\pi}{a^3}} e^{b^2/4a} + \frac{b^2}{4 a^2} \sqrt{\frac{\pi}{a}} e^{b^2/4a} \; .\]

Happily, the right-hand side gives the same:

\[ \frac{\partial^2}{\partial b^2} I_0 = \frac{\partial}{\partial b} \frac{b}{4a} \sqrt{\frac{\pi}{a}} e^{b^2/4a} = \frac{1}{2} \sqrt{\frac{\pi}{a^3}} e^{b^2/4a} + \frac{b^2}{4 a^2} \sqrt{\frac{\pi}{a}} e^{b^2/4a} \; .\]

In closing, a note on strategy. It seems that even in the case where $$b \neq 0$$ it is easier to differentiate with respect to $$a$$ for even powers.

Most students find the study of the motion of a rigid body a difficult undertaking. There are many reasons for this and it is hard to rank which one is the greatest stumbling block. Regardless of its rank, fleshing out the relationship between the angular velocity in the body-frame and the angular in the inertial frame is certainly near the top.

First off, the very concept of angular velocity being represented by a vector is foreign concept and difficult to reconcile with the idea that finite rotations are represented by matrices that do not commute. Second, we usually want to express the angular velocity in terms of time derivatives of certain Euler angles and each of these is most naturally expressed in a frame different from all the others.

This short column is intended to give a step-by-step calculation that addresses some of the difficulties encountered in this second point. While the results are not new, the presentation is fairly novel in the pedagogy, which is clean and straightforward.

For sake of concreteness, this presentation will focus on the traditional Euler 3-1-3 sequence used in celestial mechanics to orient the plane of an orbit. This selection has two advantages: it is commonly found in textbooks (e.g. Goldstein) and it involves only two transformation primitives. Its one drawback is that modern attitude control typically uses a different sequence. This latter point may be address in a companion entry at a future date.

We will use the transformation primitives

\[ {\mathcal T}_x(\alpha) = \left[ \begin{array}{ccc} 1 & 0 & 0 \\ 0 & \cos \alpha & \sin \alpha \\ 0 & -\sin \alpha & \cos \alpha \end{array}\right] \; \]

and

\[ {\mathcal T}_z(\alpha) = \left[ \begin{array}{ccc} \cos \alpha & \sin \alpha & 0 \\ -\sin \alpha & \cos \alpha & 0 \\ 0 & 0 & 1 \end{array}\right] \; . \]

The classical Euler 3-1-3 sequence gives a transformation matrix from inertial to body frame expressed as

\[ {\mathcal A} = {\mathcal T}_z(\psi) {\mathcal T}_x(\theta) {\mathcal T}_z(\phi) \; . \]

Before finding the angular velocity in the body frame, an explicit form of the attitude, $${\mathcal A}$$, is needed as it serves two purposes. It is the operator that maps body rates in the inertial frame to the body frame. It is also the device by which we get the intermediate frame orientations during the transformation. These intermediate frames provide the most straightforward route for expressing the Euler angle rates as inertial angular velocity vectors.

Start with a a known inertial frame here denoted as $$\left\{\hat x_{GCI}, \hat y_{GCI}, \hat z_{GCI} \right\}$$

and begin by applying the operator $${\mathcal T}_z(\phi)$$ corresponding to a tranformation into a frame rotated by an angle $$\phi$$ about $$\hat z_{GCI}$$.

The intermediate transformation matrix

\[ {\mathcal A}^{(1)} = \left[ \begin{array}{ccc} \cos \phi & \sin \phi & 0 \\ -\sin \phi & \cos \phi & 0 \\ 0 & 0 & 1 \end{array}\right] \; \]

allows us to read off the basis vectors in the intermediate frame from the rows of $${\mathcal A}^{(1)}$$:

\[ \hat x^{(1)} = \left[ \begin{array}{c} \cos \phi \\ \sin \phi \\ 0 \end{array} \right]_{GCI} \; ,\]

\[ \hat y^{(1)} = \left[ \begin{array}{c} -\sin \phi \\ \cos \phi \\ 0 \end{array} \right]_{GCI} \; ,\]

and

\[ \hat z^{(1)} = \left[ \begin{array}{c} 0 \\ 0 \\ 1 \end{array} \right]_{GCI} \; .\]

The subscript $$GCI$$ reminds us that these vectors are expressed in the inertial frame even though they point along the axes of the intermediate frame.

Next apply the operator $${\mathcal T}_x(\theta)$$ corresponding to a tranformation into a frame rotated by an angle $$\theta$$ about $$\hat x ^{(1)}$$.

The corresponding intermediate transformation matrix

\[ {\mathcal A}^{(2)} = \left[ \begin{array}{ccc} \cos \phi & \sin \phi & 0 \\ – \sin \phi \cos \theta & \cos \phi \cos \theta & \sin \theta \\ \sin \phi \sin \theta & -\cos \phi \sin \theta & \cos \theta \end{array}\right] \; \]

and the basis vectors in this frame are

\[ \hat x^{(2)} = \left[ \begin{array}{c} \cos \phi \\ \sin \phi \\ 0 \end{array} \right]_{GCI} \; ,\]

\[ \hat y^{(2)} = \left[ \begin{array}{c} -\sin \phi \cos \theta \\ \cos \phi \cos \theta \\ \sin \theta \end{array} \right]_{GCI} \; ,\]

and

\[ \hat z^{(2)} = \left[ \begin{array}{c} \sin \phi \cos \theta \\ -\cos \phi \sin \theta \\ \cos \theta \end{array} \right]_{GCI} \; .\]

The last step is to appluy $${\mathcal T}_z(\psi)$$ corresponding to a transformation into a frame rotated by an angle $$\psi$$ about $$\hat z^{(2)}$$.

The final transformation matrix is one of the ingredients we were seeking as is given by

\[ {\mathcal A} = \left[ \begin{array}{ccc} \cos \phi \cos \psi – \sin \phi \sin \psi \cos \theta & \cos \phi \sin \psi \cos \theta + \sin \phi \cos \psi & \sin \psi \sin \theta \\ -\sin \phi \cos \psi \cos \theta – \cos \phi \sin \psi & \cos \phi \cos \psi \cos \theta – \sin \phi \sin \psi & \cos \psi \sin \theta \\ \sin \phi \sin \theta & -\cos \phi \sin \theta & \cos \theta \end{array}\right] \; . \]

Each of the Euler angles has an associated rate $$\left\{\dot \phi, \dot \theta, \dot \psi\right\}$$. The angular velocity in the body frame is related to these rates through

\[ \vec \omega = \vec \omega_{\phi} + \vec \omega_{\theta} + \vec \omega_{\psi} = \dot \phi \hat z_{GCI} + \dot \theta \hat x^{(1)} + \dot \psi \hat z^{(2)} \; .\]

The utility of pulling out the intermediate basis vectors should now be clear.

Carrying out the multiplication gives

\[ \vec \omega = \left[ \begin{array}{c} \dot \theta \cos \phi + \dot \psi \sin \phi \sin \theta \\ \dot \theta \sin \phi – \dot \psi \cos \phi \sin \theta \\ \dot \phi + \psi \cos \theta \end{array} \right]_{GCI} \; ,\]

which is the accepted answer (e.g. problem 19 in Chapter 4 of Goldstein).

The final step is to express this vector’s components in the body frame. This is done by multiplying the column array of inertial components by $${\mathcal A}$$ to get

\[ \vec \omega = \left[ \begin{array}{c} \dot \phi \sin \psi \sin \theta + \dot \theta \cos \psi \\ \dot \phi \cos \psi \sin \theta – \dot \theta \sin \psi \\ \dot \phi \cos \theta + \dot \psi \end{array} \right]_{Body} \; ,\]

which, again, can be favorably compared with Section 4-9 of Golstein (page 179 in the second edition).

The method outlined above is the most straightforward and pedagogically sound one I know. Admittedly, the multiplications and trigonometric simplifications are tedious but with a modern computer algebra system, such as wxMaxima, it can be readily done.

The world of particle physics abounds with a lot of convenient but specialized jargon that takes some getting used too. This is particularly true in the context of natural units where the speed of light $$c$$ is a fundamnetal yark stick for measuring velocity and mass/energy and where electric charge is measured in units of the charge of an electron. The use of natural units traces back to the roots of the discipline in the studies into electromagnetism. Those practiced in the art specfy energies and masses in units of electron-volts ($$eV$$) and electron-volts per $$c^2$$.

To put things in perspective, note that room temperature $$T \approx 300 k$$ corresponds to an energy content of about $$1/40 eV$$. In contrast, the energy content associated with the mass of an electron is approximately $$511,000 eV$$ or $$511 KeV$$. The key for determining when particle motion can be treated classically versus relativistically is the ratio of kinetic energy to the rest mass energy. At the low energies representative of routine atomic processes ($$< 100 eV$$), the kinetic energy is very small compared with the rest energy and classical methods work well. On the other side are ultra-relativistic energies (such is seen in the LHC) where the kinetic energy is so much larger than the rest energy that the latter can be ignored. The grey area is that stretch inbetween low and ultra-high energy where the ratio ranges from 0.1 to 10 or so. The aim of this note is to look a little at this inbetween zone where the classical description breaks down and gives way to the relativistic one. To set the notation for the relativistic descripion, note that the interval is \[ c^2 dt^2 - dx^2 - dy^2 - dz^2 = c^2 d\tau^2 \] and that coordinate velocities are denoted as $$v^x \equiv \frac{dx}{dt} $$ and so on. Dividing out $$c dt$$ from all terms in the interval gives \[ c^2 dt^2 \left[ 1 - \frac{\vec v^2}{c^2} \right] = c^2 d\tau^2 \; ,\] from which the traditional, relativistic $$\gamma$$ is defined as: \[ \left(\frac{dt}{d\tau}\right)^2 = \frac{1}{\left[1-\vec v^2/c^2\right]} \equiv \gamma^2 \; .\] If a differential movement in spacetime is denoted by \[ dq^{\mu} = \left[\begin{array}{cccc} c dt & dx & dy & dz \end{array} \right] \] and the Minkowski metric \[ \eta_{\mu\nu} = diag(1,-1,-1,-1) \] is used, then the four-velocity is given by \[ u^{\mu} \equiv \frac{d q^{\mu}}{d \tau} = \left[\begin{array}{cccc} c \frac{dt}{d\tau} & \frac{dx}{d\tau} & \frac{dy}{d\tau} & \frac{dz}{d\tau} \end{array} \right] \; .\] The chain rule proves particularly helpful, allowing the coordinate time $$t$$ to be used in place of the proper time $$\tau$$:0 \[ u^{\mu} = \frac{dt}{d\tau}\left[\begin{array}{cccc} c & \frac{dx}{d\tau} \frac{d\tau}{dt}& \frac{dy}{d\tau}\frac{d\tau}{dt} & \frac{dz}{d\tau}\frac{d\tau}{dt} \end{array} \right] \; .\] Factoring out the $$\gamma$$ gives a compact form for the four-velocity of \[ u^{\mu} = \gamma \left[\begin{array}{cc}c & \vec v \end{array} \right] \; .\] The four-momentum is defined as \[ p^{\mu} = m u^{\mu} = m \gamma \left[\begin{array}{cc}c & \vec v \end{array} \right] \; . \] The normalization of the four-momentum is of particular physical interest. It is defined as: \[ p^{\mu} p_{\mu} = m^2 \gamma^2 \left[\begin{array}{cc}c & \vec v \end{array} \right] \left[\begin{array}{cccc}1 & 0 & 0 & 0 \\ 0 & -1 & 0 & 0 \\ 0 & 0 & -1 & 0 \\ 0 & 0 & 0 & -1 \end{array} \right] \left[\begin{array}{c}c \\ \vec v \end{array} \right] \; ,\] which simplifies to \[ p^{\mu}p_{\mu} = m^2 \gamma^2 (c^2 - \vec v^2) = m^2 \gamma^2 c^2 \left[1-\vec v^2/c^2\right] = m^2 c^2 \; . \] The zeroth component, \[ p^0 = m c \gamma = \frac{m c}{\sqrt{1-\frac{\vec v^2}{c^2} }} \; ,\] is related to the total energy by the equation \[ E = c p^0 = mc^2 \gamma = \frac{RE}{\sqrt{1 - \frac{\vec v^2}{c^2}}} \; .\] This last relation is cornerstone for looking at the cross-over from classical to relativistic descriptions. The approach will be to compute the coordinate velocity both ways and to compare the relative error made in the classical expression. Two different scenarios will be explored. The first is the usual one where we know the kinetic energy of the particle, having accelerated it through fixed potential difference. The second, and far less common scenario, is where we know the relativistic momentum of the particle after a collision. Scenario 1: Measure the kinetic energy By its thermodynamic definition, energy is additive and the total energy of a particle is given by the sum \[ E = K + RE \] of its kinetic energy $$K$$ and rest energy $$RE = m c^2$$. Recasting the definition of energy given above in terms rest energy gives \[ E = \frac{RE}{\sqrt{1 - \frac{\vec v^2}{c^2}}} \; .\] Equating the two expressions gives \[ \frac{RE}{\sqrt{1 - \frac{\vec v^2}{c^2}}} = K + Re \; ,\] which can be simplified to \[ K = (\gamma - 1) RE \; .\] Next define the ratio of kinetic and rest energies as \[ b = \frac{K}{RE} \; .\] Using this definition leads to the identity \[ b + 1 = \gamma \] and the subsequent expression for the coordinate velocity \[ v = c \sqrt{ 1 + \frac{1}{(1+b)^2} } \; \] For a $$K = 30 \; KeV$$ electron, with rest energy $$RE = 510.839 \; KeV$$, the resulting coordinate velocity (really speed) is \[ v = 98458.806 \; km/s \; .\] The classical method, expressed via the usual relation $$K = 1/2 m v^2$$, yields \[ v_{classical} = \sqrt{\frac{2 K}{m_e}} \] or, numerically, \[ v_{classical} = 102743.45 \; km/s \; ,\] corresponding to a relative error of $$\epsilon_1 = 0.0435171$$. While not a gross error, using the classical expression results in a non-negligible difference. Scenario 2: Measure the relativistic momentum If instead of knowing the kinetic energy, suppose we know the spatial portion of the relativistic four-momentum $${\mathcal P}^2 = \vec p^2 c^2$$. The normaliztion of the four-momentum means that \[ E^2 = m^2 c^2 + p^2 c^2 \equiv RE^2 + {\mathcal P}^2 \; .\] Defining the ratio \[ d = \frac{{\mathcal P}}{RE} \] and factoring and taking the square root gives \[ E = RE \sqrt{1 + d^2} \; .\] Combining with the formula for $$E$$ in terms of $$\gamma$$ allows us to determine \[ v = c\sqrt{1-\frac{1}{1+d^2}} = c \sqrt{\frac{d^2}{1+d^2}} \; .\] If we suppose that the energy associated with $${\mathcal P} = 30 KeV$$ then the resulting speed is \[ v = 17575.59 \; km/s \] whereas the classical expression is \[ v_{classical} = 17598.29 \; km/s \; .\] The relative error is $$\epsilon_2 = 0.00192$$; substantially smaller than the case reflecting the fact that the speed is so much smaller (almost a factor of 6 smaller). It is interesting to note that, in both cases, the classical expression overestimates the speed. This is consistent with the notion that the classical theory has no fundamental speed limit like relativity but it this is only a conjecture at this point.