Quasi-linear Theory: Buneman Instability¶
J. Pavan, P. H. Yoon and T. Umeda, 2011, Physics of Plasma
We reproduce results shown in the reference paper in this post. The purpose of the paper is to calculate the early stage evolution of a velocity distribution function due to Buneman instability in the framework of the quasi-linear theory.
We solve the following equations,
Here, the first equation is the electrostatic dispersion relation, second the quasi-linear diffusion equation, and thrid the equation of the spectral density of the electric field.
The coefficient, $D(v,t)$, is a diffusion coefficient which is determined as follows,
$$
D(v,t)=\frac{2}{\epsilon_0}\left(\frac{e}{m_e}\right)^2\int_{-\infty}^{\infty}\frac{\mathcal{E}(k,t)\gamma(k,t)}{\{\omega_r(k,t)-kv\}^2+\gamma^2(k,t)}dk.
$$
Note that the distribution function, $f_e(v,t)$, is spatially averaged.
The quasi-linear diffusion equation can be discretized as follows, $$ \frac{\left<f_e\right>_i^{n+1}-\left<f_e\right>_i^n}{\Delta t}=\frac{F_{i+1/2}^{n+1}-F_{i-1/2}^{n+1}}{\Delta v}; $$
$$ F_{i+1/2}^{n+1}=D_{i+1/2}^n\frac{\left<f_e\right>_{i+1}^{n+1}-\left<f_e\right>_i^{n+1}}{\Delta v}. $$Then, we solve the equation with an implicit manner.
The spectral density equation is approximated as $$ \mathcal{E}(k,t)\propto e^{2\gamma(k,t)t}. $$ Here, it is assumed that the growth rate is much slower than the time scale of the spectral density.
Initial Condition¶
The initial distribution function of electron is
$$
f_e(v)=\frac{e^{(v-v_d)^2/2v_{te}^2}}{\sqrt{2\pi}v_{te}},
$$
we set $v_{te}=1$, and $v_d=4v_{te}$.
The ion to electron mass ratio is $m_i/m_e=1600$.
The wavenumber range is $0.01\le kv_{te}/\omega_{pe}\le0.35$ with $500$ grid points, and the velocity space range is $-6\le v/v_{te}\le 10$ with $500$ grid points.
Dispersion Relation¶
%config InlineBackend.figure_format='retina'
import plasmapy.dispersion.dispersionfunction as pdd
import matplotlib.pyplot as plt
import numpy as np
from scipy import optimize
def disp(w,k):
Zi=pdd.plasma_dispersion_func_deriv(w/(np.sqrt(2)*k*vti))
Di=-wpi**2/(2*k**2*vti**2)*Zi
wr=np.real(w)
wi=np.imag(w)
zetar=(wr-k*vv)/(k*eps)
zetai=wi/(k*eps)
zeta=zetar+1j*zetai
Ze=pdd.plasma_dispersion_func_deriv(zeta)
intfeZ=np.trapz(fe*Ze,vv,dv)
De=-wpe**2/(k**2*eps**2)*intfeZ
return 1+Di+De
def freq(ww):
sol=optimize.newton(lambda x: disp(x,kk[0]), x0=0)
#sol=mm.findroot(lambda x: disp(x,kk[0]),x0=0.+1e-8j, solver='muller')
for ik in range(len(kk)):
sol=optimize.newton(lambda x: disp(x,kk[ik]),x0=sol)
#sol=optimize.root(lambda x: disp(x,kk[ik]), x0=0.9*sol,method='krylov')
#sol=mm.findroot(lambda x:disp(x,kk[ik]),x0=0.9*sol, solver='muller')
ww[ik]=sol
wr=np.real(ww)
wi=np.imag(ww)
wi[wi<0]=0
ww=wr+1j*wi
return ww
Diffusion Coefficient¶
def diffc(E,ww):
wr=np.real(ww)
wi=np.imag(ww)
D=np.zeros(len(vv))
for iv in range(len(vv)):
c1 =E*wi/((wr-kk*vv[iv])**2+wi**2)
D[iv]=2*np.trapz(c1,kk)
return D
Quasi-linear Diffusion Equation¶
def diffeq(fe,D):
wacc=1.88
cc=dt/dv**2
itermax=5000
fe_old=np.copy(fe)
Dh=np.zeros(len(vv))
for iv in range(len(vv)-1):
Dh[iv]=0.5*(D[iv+1]+D[iv])
for iter in range(itermax):
resid=0
for iv in range(1,len(vv)-1):
tmp=np.copy(fe[iv])
fe[iv]=fe_old[iv]+cc*(Dh[iv]*fe[iv+1]+Dh[iv-1]*fe[iv-1])
fe[iv]=fe[iv]/(1+cc*Dh[iv]+cc*Dh[iv-1])
fe[iv]=(1-wacc)*tmp+wacc*fe[iv]
resid=resid+np.abs(fe[iv]-tmp)
### boundary condition ###
fe[0]=fe[1]
fe[len(vv)-1]=fe[len(vv)-2]
if resid <= 1e-5:
break
#print(iter,resid)
return fe
Spectral Density Equation¶
def spectralden(E,ww):
wi=np.imag(ww)
return E*np.exp(2*wi*dt)
Main Code¶
### parameter ###
nk =500
nv =500
nt =600
vte=1
wpe=1
rm =1600
wpi=wpe/np.sqrt(rm)
vti=vte/np.sqrt(rm)
vd =4
dt =1
isav=150
### velocity space range ###
vv=np.linspace(-6,10,nv)
dv=(vv[1]-vv[0])
eps=dv#/np.sqrt(np.pi)
### initial distribution function ###
fe=np.exp(-(vv-vd)**2/(2*vte**2))/(np.sqrt(2*np.pi)*vte)
### wavenumber and frequency range
kk=np.linspace(1e-2,0.35,nk)
ww=np.zeros(nk,dtype='complex')
### initial spectral density ###
E=1e-5*np.ones(nk)
### prepare output data ###
fehst=np.zeros((1+nt//isav,nv))
wwhst=np.zeros((1+nt//isav,nk),dtype='complex')
Ehst =np.zeros((1+nt//isav,nk))
fehst[0,:]=fe
wwhst[0,:]=ww
Ehst[0,:] =E
### iterations start ###
for it in range(nt+1):
ww=freq(ww)
D =diffc(E,ww)
fe=diffeq(fe,D)
E =spectralden(E,ww)
if(it%isav==0):
print(it)
fehst[it//isav,:]=fe
wwhst[it//isav,:]=ww
Ehst[it//isav ,:] =E
Plot¶
fig, ax = plt.subplots(nrows=2,ncols=2,figsize=(16, 8))
for i in range(5):
ax[0,0].plot(kk,np.imag(wwhst[i,:]),label=r'$\omega_{pe}T=%d$' %(i*150))
ax[1,0].plot(kk,np.real(wwhst[i,:]))
ax[0,1].plot(vv,fehst[i,:])
ax[1,1].plot(kk,Ehst[i,:])
ax[0,0].legend()
ax[0,0].set_xlim( 0,0.35)
ax[1,0].set_xlim( 0,0.35)
ax[1,1].set_xlim( 0,0.35)
ax[0,1].set_xlim(-6,10 )
ax[0,0].set_ylim(0,0.05)
ax[1,0].set_ylim(0,0.07)
ax[0,1].set_ylim(0,0.4)
ax[1,1].set_ylim(0,6)
ax[0,0].set_xlabel(r'$kv_{te}/\omega_{pe}$')
ax[1,0].set_xlabel(r'$kv_{te}/\omega_{pe}$')
ax[0,1].set_xlabel(r'$v/v_{te}$')
ax[1,1].set_xlabel(r'$kv_{te}/\omega_{pe}$')
ax[0,0].set_ylabel(r'$\gamma/\omega_{pe}$')
ax[1,0].set_ylabel(r'$\omega_r/\omega_{pe}$')
ax[0,1].set_ylabel(r'$f_e$')
ax[1,1].set_ylabel(r'$\mathcal{E}$')
plt.show()
Note: the overall evolution of the distribution function and the spectral density is consistent with results in the reference paper. However, here we arbitrary ignored negative growth rate by making it zero. This treatment may not be appropriate, and have to improve the root finding technique.
!jupyter nbconvert --to html --template tmp.tpl QL_Buneman_instability.ipynb