Created by Rui Calado and Jorge Vieira, 2018
In this notebook, we are going to study the dispersion relation for electron plasma waves.
Electron plasma waves are longitudinal waves that may propagate in unmagnetized plasmas. To derive the dispersion relation for such waves let us start by considering the following setup:
We start by writing the continuity and momentum equations for the electron and ion species:
$$\large \left\{\begin{array}{lcr} \frac{\partial n_{e,i}}{\partial t}+\nabla\cdot(n_{e,i}\mathbf{v}_{e,i})=0 \\ \frac{\partial \mathbf{v}_{e,i}}{\partial t}=\mp \frac{e}{m_{e,i}}\mathbf{E}\\ \end{array}\right. . $$Then we consider Poisson's equation:
$$\epsilon_0\nabla\cdot\mathbf{E}=e(n_i-n_e).$$Applying a time derivative twice,
$$\epsilon_0\nabla\cdot\left(\frac{\partial^2 \mathbf{E}}{\partial t^2}\right)=e\left(\frac{\partial^2n_i}{\partial t^2}-\frac{\partial^2n_e}{\partial t^2}\right).$$Using the continuity and momentum equations we get:
$$\frac{\partial^2 \mathbf{E}}{\partial t^2}-\frac{e^2n_0}{\epsilon_0}\left(\frac{1}{m_i}+\frac{1}{m_e}\right)\mathbf{E}=0.$$This is the equation for a harmonic oscillator. Neglecting the $1/m_i$ term since $m_e\ll m_i$, our oscillation frequency is the electron plasma frequency:
$$\omega=\frac{e^2n_0}{\epsilon_0 m_e}\equiv \omega_{p}.$$Now, what happens if we consider a warm plasma instead? Neglecting ion motion, the first step is to add a pressure term to the electron momentum equation:
$$\frac{\partial \mathbf{v}_{e}}{\partial t}=- \frac{e}{m_{e}}\mathbf{E}-\frac{\gamma k_BT_e}{m_en_0}\nabla n_1.$$We also note that Poisson's equation now takes the form:
$$\nabla\cdot\mathbf{E}=-\frac{e}{\epsilon_0}n_1.$$Taking the divergence of the momentum equation, we get:
$$\nabla\cdot\left( \frac{\partial \mathbf{v}}{\partial t} \right)=-\frac{e}{m_e}\nabla\cdot\mathbf{E}-\frac{\gamma kT_e}{m_en_0}\nabla\cdot (\nabla n_e).$$Using the unchanged continuity equation and Poisson's equation:
$$\frac{\partial^2n_1}{\partial t^2}+\omega_{p}^2n_1-\frac{\gamma k_BT_e}{m_e}\nabla\cdot(\nabla n_1)=0.$$Considering the high frequency regime, there will be no heat losses in our time scale, and so we will take the adiabatic coefficient $\gamma=3$ for 1D longitudinal oscillations. Additionally, we use the definition $v_{th}^2=k_BT_e/m_e$ to write:
$$\frac{\partial^2n_1}{\partial t^2}+\omega_{p}^2n_1-3v_{th}^2\nabla\cdot(\nabla n_1)=0.$$The final step consists of considering sinusoidal waves such that $n_1=\text{n}_1\exp^{i(\mathbf{k}\cdot\mathbf{r}-\omega t)}$ and then Fourier analyzing the equation $\left(\nabla=i\mathbf{k},\ \frac{\partial}{\partial t}=-i\omega \right)$, which results in the dispersion relation:
$$\omega^2=\omega_{p}^2+3v_{th}^2k^2.$$
In [12]:
import em1ds as zpic
#v_the = 0.001
v_the = 0.02
#v_the = 0.20
electrons = zpic.Species( "electrons", -1.0, ppc = 64, uth=[v_the,v_the,v_the])
sim = zpic.Simulation( nx = 500, box = 50.0, dt = 0.0999/2, species = electrons )
sim.filter_set("sharp", ck = 0.99)
#sim.filter_set("gaussian", ck = 50.0)
We run the simulation up to a fixed number of iterations, controlled by the variable niter
, storing the value of the EM field $E_z$ at every timestep so we can analyze them later:
In [13]:
import numpy as np
niter = 4000
Ex_t = np.zeros((niter,sim.nx))
Ez_t = np.zeros((niter,sim.nx))
tmax = niter * sim.dt
print("\nRunning simulation up to t = {:g} ...".format(tmax))
while sim.t <= tmax:
print('n = {:d}, t = {:g}'.format(sim.n,sim.t), end = '\r')
Ex_t[sim.n,:] = sim.emf.Ex
Ez_t[sim.n,:] = sim.emf.Ez
sim.iter()
print("\nDone.")
In [14]:
import matplotlib.pyplot as plt
iter = sim.n//2
plt.plot(np.linspace(0, sim.box, num = sim.nx),Ex_t[iter,:], label = "$E_x$")
plt.plot(np.linspace(0, sim.box, num = sim.nx),Ez_t[iter,:], label = "$E_z$")
plt.grid(True)
plt.xlabel("$x_1$ [$c/\omega_n$]")
plt.ylabel("$E$ field []")
plt.title("$E_x$, $E_z$, t = {:g}".format( iter * sim.dt))
plt.legend()
plt.show()
To analyze the dispersion relation of the electrostatic plasma waves we use a 2D (Fast) Fourier transform of $E_x(x,t)$ field values that we stored during the simulation. The plot below shows the obtained power spectrum alongside the theoretical prediction.
Since the dataset is not periodic along $t$ we apply a windowing technique (Hanning) to the dataset to lower the background spectrum, and make the dispersion relation more visible.
In [15]:
import matplotlib.pyplot as plt
import matplotlib.colors as colors
# (omega,k) power spectrum
win = np.hanning(niter)
for i in range(sim.nx):
Ex_t[:,i] *= win
sp = np.abs(np.fft.fft2(Ex_t))**2
sp = np.fft.fftshift( sp )
k_max = np.pi / sim.dx
omega_max = np.pi / sim.dt
plt.imshow( sp, origin = 'lower', norm=colors.LogNorm(vmin = 1.0),
extent = ( -k_max, k_max, -omega_max, omega_max ),
aspect = 'auto', cmap = 'gray')
k = np.linspace(-k_max, k_max, num = 512)
w=np.sqrt(1 + 3 * v_the**2 * k**2)
plt.plot( k, w, label = "Electron Plasma Wave", color = 'r',ls = '-.' )
plt.ylim(0,2)
plt.xlim(0,k_max)
plt.xlabel("$k$ [$\omega_n/c$]")
plt.ylabel("$\omega$ [$\omega_n$]")
plt.title("Wave dispersion relation")
plt.legend()
plt.show()
To analyze the dispersion relation of the electrostatic plasma waves we use a 2D (Fast) Fourier transform of $E_z(x,t)$ field values that we stored during the simulation. The plot below shows the obtained power spectrum alongside the theoretical prediction.
Since the dataset is not periodic along $t$ we apply a windowing technique (Hanning) to the dataset to lower the background spectrum, and make the dispersion relation more visible.
In [8]:
import matplotlib.pyplot as plt
import matplotlib.colors as colors
# (omega,k) power spectrum
win = np.hanning(niter)
for i in range(sim.nx):
Ez_t[:,i] *= win
sp = np.abs(np.fft.fft2(Ez_t))**2
sp = np.fft.fftshift( sp )
k_max = np.pi / sim.dx
omega_max = np.pi / sim.dt
plt.imshow( sp, origin = 'lower', norm=colors.LogNorm(vmin = 1e-5, vmax = 0.01),
extent = ( -k_max, k_max, -omega_max, omega_max ),
aspect = 'auto', cmap = 'gray')
k = np.linspace(-k_max, k_max, num = 512)
w=np.sqrt(1 + k**2)
plt.plot( k, w, label = "$\omega^2 = \omega_p^2 + k^2 c^2$", color = 'r', ls = '-.' )
plt.ylim(0,k_max)
plt.xlim(0,k_max)
plt.xlabel("$k$ [$\omega_n/c$]")
plt.ylabel("$\omega$ [$\omega_n$]")
plt.title("EM-wave dispersion relation")
plt.legend()
plt.show()
In [ ]: