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

  • B0=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ϵ0tE1 \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ϵ0tE1] -\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

2E11c22t2E1(E1)=μ0tJ1 \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

2E11c22t2E1=μ0tJ1=μ0en0t[vive]=μ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

2E11c22t2E11c2[Ω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

2E11c22t2E11c2[ω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(kxωt)\vec{E} = \vec{\bar{E}} e^{i(\vec{k} \cdot \vec{x} - \omega t)} , we finally obtain

[k2+ωc2ωp2c2]E¯ei(kxω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 vgcv_g \rightarrow c as kk \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 t100t \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?

© 2018 Nextjournal GmbH