Micah P. Dombrowski / Jul 22 2018

Project 2: Electromagnetic 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/light-wave-dispersion

In this project, we are going to look at the dispersion relation for electromagnetic waves.

The dispersion relation,  ω(k)\omega(k), tells us the natural frequencies of oscillations for these waves, and the information contained in this function about the relationship between  ω\omega and  kk can be used to determine the phase and group velocities of these waves. [There will be a subsequent notebook on wave velocities]

For transverse waves have:

  •   ∇⋅E⃗=0\nabla \cdot \vec{E} = 0 -- transverse waves

  •   Te=0T_e = 0 -- cold plasma

  •   B⃗0=0\vec{B}_0 = 0 -- unmagnetized

From Maxwell's Equations we have:

−∇×E1⃗=∂∂tB1⃗ -\nabla\times\vec{E_1} = \frac{\partial}{\partial t} \vec{B_1}

and

∇×B1⃗=μ0J1⃗+μ0ϵ0∂∂tE1⃗ \nabla\times\vec{B_1} = \mu_0\vec{J_1} + \mu_0 \epsilon_0 \frac{\partial}{\partial t}\vec{E_1}

Taking the curl of the first equation and substituting into it the second equation, we get:

−∇×∇×E1⃗=∂∂t∇×B1⃗=∂∂t[μ0J1⃗+μ0ϵ0∂∂tE1⃗] -\nabla\times\nabla\times\vec{E_1} = \frac{\partial}{\partial t} \nabla\times\vec{B_1} = \frac{\partial}{\partial t}[\mu_0\vec{J_1} + \mu_0 \epsilon_0 \frac{\partial}{\partial t}\vec{E_1}]

Since  −∇×∇×A⃗=∇2A⃗−∇(∇⋅A⃗)-\nabla\times\nabla\times\vec{A} = \nabla^2\vec{A} - \nabla(\nabla \cdot\vec{A}), we have

∇2E1⃗−1c2∂2∂t2E1⃗−∇(∇⋅E1⃗)=μ0∂∂tJ1⃗ \nabla^2\vec{E_1} - \frac{1}{c^2}\frac{\partial^2}{\partial t^2}\vec{E_1} - \nabla(\nabla \cdot\vec{E_1}) = \mu_0\frac{\partial}{\partial t}\vec{J_1}

For transverse waves,  ∇⋅E1⃗=0 \nabla \cdot\vec{E_1} = 0 , so

∇2E1⃗−1c2∂2∂t2E1⃗=μ0∂∂tJ1⃗=μ0en0∂∂t[vi⃗−ve⃗]=μ0e[eME1⃗+emE1⃗]ϵ0ϵ0=μ0ϵ0[e2n0ϵ0M+e2n0ϵ0m]E1⃗\begin{aligned} \nabla^2\vec{E_1} - \frac{1}{c^2}\frac{\partial^2}{\partial t^2}\vec{E_1} & = \mu_0 \frac{\partial}{\partial t}\vec{J_1} \\ & = \mu_0 e n_0 \frac{\partial}{\partial t}[\vec{v_i}-\vec{v_e}] \\ & = \mu_0 e \Big[\frac{e}{M}\vec{E_1} + \frac{e}{m}\vec{E_1}\Big]\frac{\epsilon_0}{\epsilon_0} \\ & = \mu_0 \epsilon_0 \Big[\frac{e^2 n_0}{\epsilon_0 M} + \frac{e^2 n_0}{\epsilon_0 m}\Big]\vec{E_1} \end{aligned}

Where in the third line we used Euler's equations. Plugging in our definitions for  Ωp\Omega_p and  ωp\omega_p and moving everything to the left hand side, we finally have

∇2E1⃗−1c2∂2∂t2E1⃗−1c2[Ωp2+ωp2]E⃗=0. \nabla^2\vec{E_1} - \frac{1}{c^2}\frac{\partial^2}{\partial t^2}\vec{E_1} - \frac{1}{c^2} [\Omega_p^2 + \omega_p^2] \vec{E} = 0 .

Note as in longitudinal waves in cold unmagnetized plasmas the term  [Ωp2+ωp2][\Omega_p^2 + \omega_p^2] appears. As before, we note that this is a high frequency wave and hence approximate the ions as fixed due to their large mass, hence  [Ωp2+ωp2]≃ωp2[\Omega_p^2 + \omega_p^2] \simeq \omega_p^2, and we write the above equation as

∇2E1⃗−1c2∂2∂t2E1⃗−1c2[ωp2]E⃗=0 \nabla^2\vec{E_1} - \frac{1}{c^2}\frac{\partial^2}{\partial t^2}\vec{E_1} - \frac{1}{c^2} [\omega_p^2] \vec{E} = 0

Next, assuming  E⃗=E¯⃗ei(k⃗⋅x⃗−ωt)\vec{E} = \vec{\bar{E}} e^{i(\vec{k} \cdot \vec{x} - \omega t)} , we finally obtain

[−k2+ωc2−ωp2c2]E¯⃗ei(k⃗⋅x⃗−ωt)=0 \Big[-k^2 + \frac{\omega}{c^2} - \frac{\omega_p^2}{c^2}\Big] \vec{\bar{E}} e^{i(\vec{k} \cdot \vec{x} - \omega t)} = 0
⇒ω2=ωp2+c2k2,\Rightarrow \omega^2 = \omega_p^2 + c^2k^2,

the dispersion relation for an electromagnetic wave in an unmagentized plasma!

It is plotted below.

# Plotting w(k)
import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline

N = 5
k = np.linspace(-N,N,N*20)
w_p = 1
c = 1
w = np.sqrt(w_p**2 + c**2 * k**2)
cline = k
plt.plot(k,w, label='$\omega$(k)')
plt.plot(k,cline, label='slope = c')
plt.xlabel('k [$c/\omega_p$]')
plt.ylabel('$\omega$ (in units of $\omega_p$)')
plt.xlim(-N,N)
plt.ylim(0,N+1)
plt.grid(b=True, which='major', axis='both')
plt.legend(loc=0)
plt.show()

A couple of things to note about this plot:

  • recalling expressions for phase velocity,  vϕ=ω/kv_{\phi}=\omega/k, and group velocity,  vg=dωdkv_g = \frac{d\omega}{dk} we obtain
vg=c1−ωp2/ω2 v_g = c\sqrt{1-\omega_p^2/\omega^2}
vϕ=c/1−ωp2/ω2 v_{\phi} = c / \sqrt{1-\omega_p^2/\omega^2}

  • Comparing this with the plot, confirm that indeed  vg=0v_g = 0 for  k=0k = 0, and that  vg→cv_g \rightarrow c as  k→∞k \rightarrow \infty. Importantly, although  vϕ>cv_{\phi} > c, we have  vg<cv_g < c. Thus special relativity is not violated, since information can only propagate at the group velocity and not at the phase velocity.

  • Also note, if  ω<ωp\omega < \omega_p a wave cannot propogate since  kk becomes imaginary and we get an evanescent wave.

1.
Simulations with a Particle-in-Cell Code

In this project you simulate plasmas with exactly the same conditions as in Project 1a.

Each plasma electron is initialized with positions (only in  xx or what we call  x1x_1) such that the density is uniform. The ions are initialized at the same positions but they have an infinite mass. Each electron is also initialized with velocities ( v1v_1,  v2v_2,  v3v_3) or momentum ( mv1mv_1,  mv2mv_2,  mv3mv_3) from a Maxwellian in each direction. The particles then begin to move in the self-consistent fields that their current and charge density produce.

  • The length of the plasmas is 50  c/ωpc/\omega_p 
  • The simulation will run for a time 400  1/ωp1/\omega_p.
  • The simulation uses 50,000 particles.

You will be looking at plots of the electric field in the  x3x_3 direction,  E3E_3. In Project 1a you plotted  E1E_1.

1.1.
The following lines must always be executed before running anything else.

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

import osiris
%matplotlib inline

2.
Run a case in which  vth1=vth2=vth3=0.02cvth_1=vth_2=vth_3=0.02 c.

# vth/c = 0.02
dirname = 'v02'
osiris.runosiris(rundir=dirname,inputfile='v02.txt')

After the simulation is finished, plot  E3(x1)E_3(x_1) at  t≈100t \approx 100 (run the next cell).

  • Do you see any evidence of an electromagnetic wave?
  • Does the plot make sense?
dirname = 'v02'
osiris.field(rundir=dirname, dataset='e3', time=100)

Next, plot  E3(t)E_3(t) at  x1=5c/ωpx_1=5 c/\omega_p (i.e., at cell=100).

  • Do you see any evidence of an electromagnetic wave?
  • Does the plot make sense?
dirname = 'v02'
osiris.field(rundir=dirname, dataset='e3', space=100)

Next, in the following two cells, we are going to plot  ω(k)\omega(k). This is generated by taking  E3(x1,t)E_3(x_1,t) and Fourier analyzing in both position and time.

  •   ω(k)\omega(k) with wavenumber in units of [k] =  ωpe/c\omega_{pe}/c:
dirname = 'v02'
osiris.plot_wk(rundir=dirname, wlim=[0,10], klim=[0,10], vth = 0.02, vmin=-7, vmax=2, plot_or=3) 
osiris.plot_wk(rundir=dirname, wlim=[0,10], klim=[0,10], vth = 0.02, vmin=-7, vmax=2, plot_or=3, show_theory=True)
  •   ω(k)\omega(k) with wavenumber in units of [k] =  λDe\lambda_{De}:
dirname = 'v02'
osiris.plot_wk(rundir=dirname, wlim=[0,10], klim=[0,0.4], vth = 0.02, vmin=-7, vmax=2, show_theory=True, 
               debye=True, plot_or=3)

I would like you to think about units:

  • What do you think are the natural units for  ω\omega?
  • What about k?

We are plotting them in what are called normalized units (not inverse time or inverse distance). We use two choices ( ω/ωp\omega/\omega_p and  kc/ωpkc/\omega_p) and ( ω/ωp\omega/\omega_p and  kvth/ωp=kλDkv_{th}/\omega_p = k\lambda_D).

We also plot the theory curve:

  • Does it make sense?
  • Why do you think it agrees better for all wavenumber than did the plot in Project 1a whichw as for Bohm-Gross waves?

3.
Run a case in which  vth1=vth2=vth3=0.05cvth_1=vth_2=vth_3=0.05 c.

# vth/c = 0.05
dirname = 'v05'
osiris.runosiris(rundir=dirname,inputfile='v05.txt')

Make  ω(k)\omega(k) plots for this case by running the cells below.

  •   ω(k)\omega(k) with wavenumber in units of [k] =  ωpe/c\omega_{pe}/c:
dirname = 'v05'
osiris.plot_wk(rundir=dirname, wlim=[0,10], klim=[0,10], vth = 0.05, vmin=-7, vmax=2, plot_or=3) 
osiris.plot_wk(rundir=dirname, wlim=[0,10], klim=[0,10], vth = 0.05, vmin=-7, vmax=2, plot_or=3, show_theory=True)
  •   ω(k)\omega(k) with wavenumber in units of [k] =  λDe\lambda_{De}:
dirname = 'v05'
osiris.plot_wk(rundir=dirname, wlim=[0,10], klim=[0,0.4], vth = 0.05, vmin=-7, vmax=2, show_theory=True, 
               plot_or=3, debye=True)

3.1.
Questions

  • Do the  ω(k)\omega(k) plots make sense?
  • For which normalized units do the plots look similar to case b?

4.
Run a case in which  vth1=vth2=vth3=0.20cvth_1=vth_2=vth_3=0.20 c.

# vth/c = 0.2
dirname = 'v20'
osiris.runosiris(rundir=dirname,inputfile='v20.txt')

Make  ω(k)\omega(k) plots for this case by running the cells below.

  •   ω(k)\omega(k) with wavenumber in units of [k] =  ωpe/c\omega_{pe}/c:
dirname = 'v20'
osiris.plot_wk(rundir=dirname, wlim=[0,10], klim=[0,10], vth = 0.20, vmin=-5, vmax=4, plot_or=3) 
osiris.plot_wk(rundir=dirname, wlim=[0,10], klim=[0,10], vth = 0.20, vmin=-5, vmax=4, plot_or=3, show_theory=True)
  •   ω(k)\omega(k) with wavenumber in units of [k] =  λDe\lambda_{De}:
dirname = 'v20'
osiris.plot_wk(rundir=dirname, wlim=[0,10], klim=[0,0.4], vth = 0.20, vmin=-5, vmax=4, show_theory=True, plot_or=3, debye=True)

4.1.
Questions

  • Do the  ω(k)\omega(k) plots make sense?

  • For which normalized units do the plots look similar to case a and b?

  • Look closely at  k=0k=0. The frequency does not agree with theory as well as for cases a and b.

  • Is it higher or lower than theory?

  • Can you think of a reason why?

  • Does the agreement with the theory depend much the temperature? Can you justify your answer using some equations?