Buneman Instability

Governing Equations

\begin{gathered} {\partial _t}{n_s} + \nabla \cdot \left( {{n_s}{{\mathbf{U}}_s}} \right) = 0; \hfill \\ {m_s}{n_s}\left( {{\partial _s}{{\mathbf{U}}_s} + {{\mathbf{U}}_s}\nabla {{\mathbf{U}}_s}} \right) = {n_s}{e_s}{\mathbf{E}}; \hfill \\ \nabla \cdot {\mathbf{E}} = 4\pi {\rho _q} \hfill \\ {\rho _q} = \sum\limits_s {{n_s}{e_s} } \hfill \\ \end{gathered}

1D System

\begin{gathered} {\partial _t}{n_s} + {\partial _x}\left( {{n_s}{U_{sx}}} \right) = 0; \hfill \\ {m_s}{n_s}\left( {{\partial _t}{U_{sx}} + {U_{sx}}{\partial _x}{U_{sx}}} \right) = {n_s}{e_s}{E_x}; \hfill \\ {\partial _x}{E_x} = 4\pi {\rho _q} \hfill \\ {\rho _q} = \sum\limits_s {{n_s}{e_s}; } \hfill \\ \end{gathered}

Linear Dispersion Relation

Let ${U_{ix}} = U_0+{\tilde U_i}$, ${U_{ex}} = {\tilde U_e}$, ${n_{i}} = n_{i0}+{\tilde n_i}$, ${n_{e}} = n_{e0}+{\tilde n_e}$, and ${E_{x}} = {\tilde E_x}$, we fourier-transform those equations and get the linear dispersion relation,

$$1 - \frac{{\omega _{pi}^2}}{{{{\left( {\omega - k{U_0}} \right)}^2}}} - \frac{{\omega _{pe}^2}}{{{{{\omega }}^2}}} = 0$$

%config InlineBackend.figure_format = 'retina'
import numpy as np
import matplotlib.pyplot as plt
import math as mt
from numba.decorators import jit
from numba import f8,u8
import scipy.fftpack as sf
from IPython import display

An Example

• $m_i/m_e=400$
• $U_0=0.5$

nk=100
rm=400
wpe2=1
wpi2=wpe2/rm
v0=0.5
k=np.linspace(0,4,nk)
w=np.zeros((4,nk),dtype=complex)
for ik in range(nk):
    coeff=[1.0,-2.0*k[ik]*v0,(k[ik]*v0)**2.0-wpi2-wpe2,2.0*k[ik]*v0*wpe2,-(k[ik]*v0)**2*wpe2]
    w[:,ik]=np.sort(np.roots(coeff))

plt.plot(k,np.real(w[0,:]),'k');plt.plot(k,np.imag(w[0,:]),'k--')
plt.plot(k,np.real(w[1,:]),'k');plt.plot(k,np.imag(w[1,:]),'k--')
plt.plot(k,np.real(w[2,:]),'k');plt.plot(k,np.imag(w[2,:]),'k--')
plt.plot(k,np.real(w[3,:]),'k');plt.plot(k,np.imag(w[3,:]),'k--')
plt.ylim(-2,2)
plt.xlabel(r'$k(v_0/\omega_{pe})$');plt.ylabel(r'$\omega/\omega_{pe}$')
plt.xlim(0.0,4.0)
plt.show()

Full-PIC Simulation

def es1(rm,vt,ppc,dx,nx,dt,nt,iphs):
    np.seterr(divide='ignore', invalid='ignore')
    nop=ppc*nx
    lx=dx*nx
    xx=dx*np.arange(nx)
    tt=dt*np.arange(nt)
    xp=np.zeros((2,nop))
    vp=np.zeros((2,nop))
    xp[0,:]=np.linspace(0,lx-lx/nop,nop)
    xp[1,:]=np.linspace(0,lx-lx/nop,nop)
    vp[0,:]=np.random.normal(0,vt,nop)
    vp[1,:]=np.random.normal(0,vt/mt.sqrt(rm),nop)+0.5
    
    #ds=np.zeros((2,nx))
    kk=2*mt.pi/lx*np.r_[np.arange(nx/2),np.arange(-nx/2,0)]
    Kl=(np.sin(0.5*kk*dx)/(0.5*dx))**2
    kl=np.sin(kk*dx)/dx
    
    exsav=np.zeros((nt//iphs,nx))
    nisav=np.zeros((nt//iphs,nx))
    nesav=np.zeros((nt//iphs,nx))
    drfsav=np.zeros((nt//iphs,2))
    tmpsav=np.zeros((nt//iphs,2))
    
    q,qomdt   =np.zeros((2,2))
    q[0]=-nx/nop;qomdt[0]=-dt 
    q[1]=-q[0];  qomdt[1]= dt/rm
    
    for it in range(nt):
        xp=xp+dt*vp
        xp[xp>lx]=xp[xp>lx]-lx
        xp[xp<0 ]=xp[xp<0 ]+lx
        
        ds=np.zeros((2,nx))
        ds=dens(dx,nx,nop,q,xp,ds)
        exfft=-1j/kk*sf.fft(np.sum(ds,0))
        #exfft=-1j*kl*sf.fft(np.sum(ds,0))/Kl
        exfft[0]=0
        ex=np.real(sf.ifft(exfft))        
        
        vp=acc(dx,nx,nop,xp,vp,ex,qomdt)
        
        if(it%iphs==0):
            exsav[it//iphs,:]=ex[:]
            nesav[it//iphs,:]=ds[0,:]
            nisav[it//iphs,:]=ds[1,:]
            drfsav[it//iphs,0] =np.average(vp[0,:])
            drfsav[it//iphs,1] =np.average(vp[1,:])
            tmpsav[it//iphs,0] =np.average((vp[0,:]-np.average(vp[0,:]))**2)
            tmpsav[it//iphs,1] =np.average((vp[1,:]-np.average(vp[1,:]))**2)
        
        #if(it%16==0):
        #    display.clear_output(True)
        #    plt.figure(figsize=(12,6))
        #    plt.subplot(211);plt.plot(xp[0,:],vp[0,:],',')
        #    plt.subplot(212);plt.plot(xp[1,:],vp[1,:],',')
        #    plt.show()

    return locals()

@jit('f8[:,:](f8,u8,u8,f8[:,:],f8[:,:],f8[:],f8)')
def acc(dx,nx,nop,xp,vp,ex,qomdt):
    for isp in range(2):
        for ip in range(nop):
            ixm=mt.floor(xp[isp,ip]/dx); ixp=ixm+1
            wxp=xp[isp,ip]/dx-ixm; wxm=1-wxp

            if ixp>nx-1: ixp=ixp-nx
            vp[isp,ip]=vp[isp,ip]+qomdt[isp]*(wxm*ex[ixm]+wxp*ex[ixp])             
    return vp

@jit('f8[:,:](f8,u8,u8,f8[:],f8[:,:],f8[:,:])')
def dens(dx,nx,nop,q,xp,ds):
    for isp in range(2):
        for ip in range(nop):
            ixm=mt.floor(xp[isp,ip]/dx); ixp=ixm+1
            wxp=xp[isp,ip]/dx-ixm; wxm=1-wxp

            if ixp>nx-1: ixp=ixp-nx

            ds[isp,ixm]=ds[isp,ixm]+wxm
            ds[isp,ixp]=ds[isp,ixp]+wxp 
        ds[isp,:]=ds[isp,:]*q[isp]
    return ds

%%time 
rm=400;vt=0.05;ppc=1000;nx=1024;nt=2048;iphs=1;dt=0.1;dx=0.1
data=es1(rm,vt,ppc,dx,nx,dt,nt,iphs) #(rm,vt,nop,dx,nx,dt,nt,iphs)
locals().update(data)

CPU times: user 3min 44s, sys: 45.6 s, total: 4min 30s
Wall time: 2min 23s

Space time evolution of $E_x$ and the dispersion relation

kmin=2*mt.pi/(dx*nx)*(-nx/2); kmax=2*mt.pi/(dx*nx)*(nx/2)
wmin=2*mt.pi/(dt*nt)*(-nt/2); wmax=2*mt.pi/(dt*nt)*(nt/2)
kx=np.linspace(kmin,kmax,nx); waxis=np.linspace(wmin,wmax,nt)

plt.figure(figsize=(10,5))
plt.subplot(121)
plt.imshow(exsav,origin='lower',aspect='auto',cmap='bwr',extent=([0,nx*dx,0,nt*dt]))
plt.colorbar()
plt.xlabel(r'$x/(v_0/\omega_{pe})$'); plt.ylabel(r'$t\omega_{pe}$')

plt.subplot(122)
fftex=sf.fftshift(sf.fft2(exsav[:,::-1]))
plt.imshow(abs(fftex),origin='lower',aspect='auto',extent=([kmin,kmax,wmin,wmax]))
plt.plot(k,np.real(w[0,:]),'w');plt.plot(k,np.imag(w[0,:]),'w--')
plt.plot(k,np.real(w[1,:]),'w');plt.plot(k,np.imag(w[1,:]),'w--')
plt.plot(k,np.real(w[2,:]),'w');plt.plot(k,np.imag(w[2,:]),'w--')
plt.plot(k,np.real(w[3,:]),'w');plt.plot(k,np.imag(w[3,:]),'w--')
plt.axis([0,4,-1,2])
plt.clim(0,1e4)
plt.xlabel(r'$k(v_0/\omega_{pe})$');plt.ylabel(r'$\omega/\omega_{pe}$')
plt.show()

plt.figure(figsize=(10,5))
plt.subplot(121)
plt.imshow(nisav,origin='lower',aspect='auto',cmap='Blues',extent=([0,nx*dx,0,nt*dt]))
plt.title('proton denisty')
plt.xlabel(r'$x/(v_0/\omega_{pe})$'); plt.ylabel(r'$t\omega_{pe}$')
plt.colorbar()

plt.subplot(122)
plt.imshow(nesav,origin='lower',aspect='auto',cmap='Reds',extent=([0,nx*dx,0,nt*dt]))
plt.title('electron denisty')
plt.xlabel(r'$x/(v_0/\omega_{pe})$'); plt.ylabel(r'$t\omega_{pe}$')
plt.colorbar()
plt.show()

Energy Power Spectrum

it=900
rho2fft=np.abs(sf.fftshift(sf.fft(exsav[it,:])))
plt.plot(kx,rho2fft,label=(r'$E_x$'))
plt.title(r'$t\omega_{pe}=%2.f$' %(it*dt*iphs))
plt.loglog()
plt.legend()
plt.xlabel(r'$k$'); plt.ylabel('power spectrum')
#plt.xlim(1,nx/3)
plt.ylim(1e-4,1e2); plt.yticks(np.logspace(-4,2,4))
plt.show()

Time Evolution

nd=nt//iphs
exhst=np.zeros(nd)
for it in range(nd):
    exhst[it]=np.average(exsav[it,:]**2)

tt=np.arange(nt//iphs)*(dt*iphs)
plt.plot(tt,exhst)
plt.xlim(0,nt*dt)
plt.ylabel(r'$E_x^2$')
plt.xlabel(r'$t\omega_{pe}$')
plt.semilogy()
plt.show()

fftex=np.abs(sf.fftshift(sf.fft(exsav,axis=1),axes=1))
idx=np.abs(kx-2).argmin()
tt=np.arange(nt//iphs)*(dt*iphs)
plt.plot(tt,abs(fftex[:,idx]),'-')
plt.plot(tt,1e-2*np.exp(0.09*tt),'--',label=(r'$\gamma\sim0.09$'))
plt.legend()
plt.xlabel(r'$t\omega_{pe}$')
plt.title(r'the mode: $k=2$')
plt.semilogy()
plt.show()

tt=np.arange(nt)*dt
plt.figure(figsize=(15,10))
plt.subplot(221);plt.plot(tt,drfsav[:,0],label='electron');plt.ylabel('average velocity');plt.legend()
plt.subplot(222);plt.plot(tt,drfsav[:,1],label='ion');plt.legend()
plt.subplot(223);plt.plot(tt,tmpsav[:,0]/vt**2);plt.xlabel(r'$t\omega_{pe}$');plt.ylabel(r'$T_e/T_{e0}$')
plt.subplot(224);plt.plot(tt,tmpsav[:,1]/(vt/np.sqrt(rm))**2);plt.xlabel(r'$t\omega_{pe}$');plt.ylabel(r'$T_i/T_{i0}$')
plt.show()