Micah P. Dombrowski / Jul 27 2018

Waves

Imported and adapted from the UCLA Plasma Simulation Group's JupyterPIC Notebooks 

Before anything else, we need to change to the stored notebook's directory, which contains input files for the PIC simulator.

cd /notebooks/velocities

All displays of prerendered graphics via Jupyter cell were replaced with direct uploads of the image files.

Closely related terms:

  • normal modes
  • natural oscillations / frequencies
  • resonance
  • collective modes
  • dispersion relations

Suppose you knew the initial values for the  E⃗\vec{E} and  B⃗\vec{B} fields in a plasma at some initial time. You might wonder what the fields are at any later time.

How do you answer this?

Consider the ridiculously simple example of a spring.

(Reminder: Hit Shift+Enter to run a cell, or select the cell and click on the "Run" button in the top menu bar)

The equations are:

mdvdt=Fm \frac{dv}{dt} = F
v=dxdtv = \frac{dx}{dt}
F=−kxF = -kx

In principle you could use these equations to solve for the motion of the spring.

Or since there are three equations and three unknowns (x, v, and F), you can combine these equations into a single equation of only one quantity  xx:

d2xdt2+kmx=0\frac{d^2x}{dt^2} + \frac{k}{m}x = 0

With  ω0=k/m\omega_0 = \sqrt{k/m}, we find the harmonic oscillator equation:

d2xdt2+ω02x=0\frac{d^2x}{dt^2} + \omega_0^2 x = 0

The solution to this equation is

x=Acosω0t+Bsinω0tx = A \cos{\omega_0 t} + B \sin{\omega_0 t}

where A and B are determined from the initial conditions for  xx and  vv, and  ω0\omega_0 is the "natural" frequency.

If all we wanted was the natural frequency for a (set of coupled) equation(s), we could substitute  x=ℜ[x¯eiωt]x = \Re[\bar{x}e^{i\omega t}] into the equation(s).

For our example ( d2xdt2+ω02x=0\frac{d^2x}{dt^2} + \omega_0^2 x = 0), this gives

(−ω2+ω02)x¯eiωt=0,(-\omega^2 + \omega_0^2) \bar{x}e^{i\omega t} = 0,
ω=±ω0\omega = \pm \omega_0

We could also instead substitute  x=x¯eiωtx = \bar{x}e^{i\omega t},  v=v¯eiωtv = \bar{v}e^{i\omega t}, and  F=F¯eiωtF = \bar{F}e^{i\omega t} into the original three equations, giving

[miωv¯=F¯]eiωt[mi\omega \bar{v} = \bar{F}]e^{i\omega t}
[v¯=iωx¯]eiωt[\bar{v} = i\omega \bar{x}]e^{i\omega t}
[F¯=−kx¯]eiωt[\bar{F} = -k \bar{x}]e^{i\omega t}

This is now just a set of algebraic equations to solve (rather than a differential equation).  x¯\bar{x},  v¯\bar{v}, and  F¯\bar{F} are complex numbers, sometimes called phasors.

Solving the three equations gives

miω(iωx¯)=−kx¯mi\omega(i\omega\bar{x}) = -k\bar{x}
(−ω2+km)x¯=0\left(-\omega^2 + \frac{k}{m}\right)\bar{x} = 0

and again we find  ω=k/m\omega = \sqrt{k/m}.

This method is useful because sometimes it is not possible to obtain a single equation to solve. We can "Fourier Analyze" each equation and then solve the algebraic equations for the Fourier "amplitudes."

Note that this only works for linear equations. (Why??)

If we add an external driving force to the equations:

F=−kx+FextF = -kx + F_{ext}

then

(d2dt2+ω02)x=Fext\left(\frac{d^2}{dt^2} + \omega_0^2\right)x = F_{ext}

If  Fext=F0eiωtF_{ext} = F_0 e^{i\omega t} then when  ω→ω0\omega \rightarrow \omega_0 the amplitude of x gets very large

x¯=F0−ω2+ω02\bar{x} = \frac{F_0}{-\omega^2 + \omega_0^2}

This is called resonance: Drive something at its resonance (at its natural frequency) and it gets very large.

1.
Phase velocity

We now consider waves that can move in space as well as oscillate in time.

Consider the following partial differential equation:

(∂2∂t2−vϕ2∂2∂x2)E⃗=0\left(\frac{\partial^2}{\partial t^2} - v_{\phi}^2\frac{\partial^2}{\partial x^2}\right)\vec{E} = 0

This wave equation pops up in several areas across physics. For example, by taking the curl of Faraday's law,

∇×[−∇×E⃗=∂∂tB⃗]\nabla \times \left[- \nabla \times \vec{E} = \frac{\partial}{\partial t}\vec{B} \right]

and combining it with the time derivative of Ampere's Law in the case of zero current (like in a vacuum),

∂∂t[∇×B⃗=1c2∂∂tE⃗]\frac{\partial}{\partial t} \left[ \nabla \times \vec{B} = \frac{1}{c^2}\frac{\partial}{\partial t}\vec{E} \right]

we get

∇×∇×E⃗+1c2∂2∂t2E⃗=0\nabla \times \nabla \times \vec{E} + \frac{1}{c^2}\frac{\partial^2}{\partial t^2}\vec{E} = 0

If it is 1D-like with only variations in x (so that  ∂/∂x≠0\partial/\partial x \neq 0 but  ∂/∂y=∂/∂z=0\partial/\partial y = \partial/\partial z = 0),

(∂2∂t2−c2∂2∂x2)E⃗=0\left(\frac{\partial^2}{\partial t^2} - c^2 \frac{\partial^2}{\partial x^2}\right)\vec{E} = 0

And  cc turns out to be the velocity of these electromagnetic waves.

Returning to the above equation with  vϕv_{\phi}, let

E=E¯ei(kx−ωt)E = \bar{E}e^{i(kx-\omega t)}

then

(−ω2+k2vϕ2)E¯ei(kx−ωt)=0\left(-\omega^2 + k^2 v_{\phi}^2\right)\bar{E}e^{i(kx-\omega t)} = 0
ω=±kvϕ \omega = \pm k v_{\phi}

  vϕ=ω/kv_{\phi} = \omega/k is called the phase velocity of the wave.

#####################
# CHANGE here to make different waves and watch the velocity change
w = 5.0
k = 1.0
#####################

# Thanks to http://jakevdp.github.io/blog/2013/05/12/embedding-matplotlib-animations/ and http://louistiao.me/posts/notebooks/embedding-matplotlib-animations-in-jupyter-notebooks/
%matplotlib inline
import numpy as np
import matplotlib.pyplot as plt
from matplotlib import animation, rc
from IPython.display import HTML

# First set up the figure, the axis, and the plot element we want to animate
fig, ax = plt.subplots(); ax.set_xlim(( 0, 10)); ax.set_ylim((-2, 2)); line, = ax.plot([], [], lw=2)

# initialization function: plot the background of each frame
def init():
    line.set_data([], []); return (line,)

# animation function. This is called sequentially
def animate(i):
    x = np.linspace(0, 10, 1000); y = np.sin(k * x - w * 0.05 * i); line.set_data(x, y); return (line,)

# call the animator. blit=True means only re-draw the parts that have changed.
anim = animation.FuncAnimation(fig, animate, init_func=init,frames=int(2*np.pi/w*40), interval=20, blit=True)
HTML(anim.to_html5_video())

The "phase" is  Φ=(kx−ωt)\Phi = (kx-\omega t). A peak of the wave corresponds to one value of the phase, a trough of the wave is another phase point, and etc.

If you track a peak of the wave (a point of constant phase) in space and time to see how fast it moves, you will be measuring the phase velocity. By tracking a peak where the phase is constant:

dΦdt=0\frac{d\Phi}{dt} = 0

which means that

kdxdt−ω=0k \frac{dx}{dt} - \omega = 0

from which we get

dxdt=ωk=vϕ\frac{dx}{dt} = \frac{\omega}{k} = v_{\phi}

2.
Group velocity

In addition to  vϕv_{\phi} there are other velocities associated with a wave.

To see this, consider two waves at slightly different frequencies and wavenumbers:  (ω1,k1)(\omega_1,k_1) and  (ω2,k2)(\omega_2,k_2):

E=Acos(k1x−ω1t)+Bcos(k2x−ω2t)E = A \cos{(k_1 x - \omega_1 t)} + B \cos{(k_2 x - \omega_2 t)}
E=AcosΦ1+BcosΦ2E = A \cos{\Phi_1} + B \cos{\Phi_2}

By utilizing the trigonometric identity that  cosAcosB=12(cos(A+B2)+cos(A−B2))\cos{A}\cos{B} = \frac{1}{2}(\cos(\frac{A + B}{2}) + \cos(\frac{A - B}{2})),

E=2Acos(Φ1+Φ22)cos(Φ1−Φ22)E = 2A \cos\left(\frac{\Phi_1 + \Phi_2}{2}\right) \cos\left(\frac{\Phi_1 - \Phi_2}{2}\right)
E=2Acos(k¯x−ω¯t)cos((Δk)x−(Δω)t2)E = 2A \cos(\bar{k}x - \bar{\omega}t) \cos\left(\frac{(\Delta k) x - (\Delta \omega) t}{2}\right)

where  k¯=(k1+k2)/2≈k1\bar{k} = (k_1+k_2)/2 \approx k_1 or  k2k_2, and  ω¯=(ω1+ω2)/2≈ω1\bar{\omega} = (\omega_1+\omega_2)/2 \approx \omega_1 or  ω2\omega_2.

#####################
# CHANGE here to make different waves,
# and DRAG the slider button to watch the waves move
w1 = 20.0
k1 = 1.0
w2 = 20.05
k2 = 1.1
# BEWARE: because the plots are being generated at discrete time steps and at a discrete range of points,
# you may find that waves "fictiously" move at different velocities if the grid or frame rate
# does not properly resolve them
#####################
%matplotlib inline
import numpy as np
import matplotlib.pyplot as plt
from ipywidgets import interact
x = np.linspace(0, 100, 1000);
def intplot(t=0):
    plt.figure(figsize=(8, 4))
    plt.plot(x,np.sin(k1 * x - w1 * t) + np.cos(k2 * x - w2 * t))
    plt.xlim(0,100)
    plt.ylim(-2,2)
    return plt
interact(intplot,t=(0,100))

The solution here consists of the product of two motions.

The quickly oscillating wave  (cos(k¯x−ω¯t))(\cos(\bar{k}x - \bar{\omega}t)) moves at the phase velocity

ω¯k¯≈ω1k1orω2k2 \frac{\bar{\omega}}{\bar{k}} \approx \frac{\omega_1}{k_1} \text{or} \frac{\omega_2}{k_2}

The envelope on top,  cos((Δkx−Δωt)/2)\cos((\Delta k x - \Delta \omega t)/2), moves at a different velocity called the group velocity

ΔωΔk≈∂ω∂k=vg \frac{\Delta\omega}{\Delta k} \approx \frac{\partial \omega}{\partial k} = v_g

  vgv_g can also be thought of as the energy transport velocity. One can "communicate" with energy, so  vg<cv_g < c in order for special relativity to be satisfied, but there is no such restraint on  vϕv_{\phi} (which can in fact be larger than  cc).

#####################
# CHANGE here to make different waves
# and NOTE the different slopes of v_phi vs v_g (can you tell which is which?)
w1 = 20.0
k1 = 1.0
w2 = 20.05
k2 = 1.1
# BEWARE: because the plots are being generated at discrete time steps and at a discrete range of points,
# you may find that waves "fictiously" move at different velocities if the grid or frame rate
# does not properly resolve them
#####################
%matplotlib inline
import numpy as np
import matplotlib.pyplot as plt
x = np.linspace(0, 100, 5000);
plotdata = np.zeros( (100,5000) )
for t in range(100):
    plotdata[t,:] = np.sin(k1 * x - w1 * t) + np.cos(k2 * x - w2 * t)
plt.figure(figsize=(8, 5))
plt.imshow(plotdata,origin='lower',aspect='auto',extent=[0,100,0,100],cmap="nipy_spectral")
plt.xlim(0,100); plt.ylim(0,100); plt.clim(-2,2); plt.colorbar(orientation='vertical'); 
plt.title('Time vs Space'); plt.xlabel('x'); plt.ylabel('t')
plt.show()

3.
Dispersion

There are two key velocities:

vϕ=ωkv_{\phi} = \frac{\omega}{k}
vg=∂ω∂kv_g = \frac{\partial \omega}{\partial k}

These velocities depend on the relationship between  ω\omega and  kk. Such a function,

ω(k)\omega(k)

is called the dispersion relation of the wave.

Pick a  kk, and then  ω\omega is defined. Pick an  ω\omega, and then  kk is defined.

We need dispersion relations to determine the behavior of waves in plasmas.

For the next several weeks, our goal will be to derive the dispersion relation  ω(k)\omega(k) for many of the basic waves which can exist in a plasma. (not all - not enough time.)