Numerical Cherenkov Instability

%%bash
cat input.f90

module input
!### definition by a user ###!
integer, parameter :: ppc=10
integer, parameter :: nx=128
integer, parameter :: ny=128
integer, parameter :: nsp=2
integer, parameter :: nt=1000
integer, parameter :: isavw=10, isavp=100
real(8), parameter :: theta=0.505d0
real(8), parameter :: c =1.d0
real(8), parameter :: dx=0.1d0
real(8), parameter :: dy=0.1d0
real(8), parameter :: dt=0.05d0
real(8), parameter :: b0=0.d0
real(8), parameter :: th=0, phi=0
real(8), parameter :: q(nsp)    =(/ -1.d0             , 1.d0               /)
real(8), parameter :: rm(nsp)   =(/  1.d0             , 1.d0               /)
real(8), parameter :: qomdt(nsp)=(/  0.5*dt*q(1)/rm(1), 0.5*dt*q(2)/rm(2)  /)
real(8), parameter :: dn(nsp)   =(/  1.d0/ppc         , 1.d0/ppc           /)
real(8), parameter :: vt(nsp)  =(/   0.1d0            , 0.1d0/sqrt(rm(2))  /)
real(8), parameter :: vx0(nsp)  =(/  100.d0           , 100.d0             /)

!### don't change if you have no reason! ###!
integer, parameter :: npmax=ppc*nx*ny
integer, dimension(nsp) :: np
integer, dimension(nx) :: iipx,iimx
integer, dimension(ny) :: iipy,iimy
integer, dimension(nx*ny) :: ijmx, ijpx, ijmy, ijpy
real(8), parameter :: pi=3.141592653589793d0
real(8), parameter :: lx=dx*nx, ly=dy*ny

end module input

Parameters

• $m_i/m_e=1$
• $ppc=10$
• $v_{th}=0.1c$
• $\Gamma=100$
• $c\Delta t/\Delta x=0.5$

%config InlineBackend.figure_format = 'retina'
from numba import jit,f8,u8
import numpy as np
import matplotlib.pyplot as plt
import math as mt
from scipy.fftpack import fft, ifft,fft2,fftshift
import time
import matplotlib as mpl
from IPython import display
plt.rcParams['font.size'] = 14
plt.rcParams['axes.linewidth'] = 1.5

nx=128
ny=128
nd=100
isav=10
dx=0.1; dy=0.1
dt=0.05
lx=nx*dx; ly=ny*dy

Snapshot in the real space

it=50
f =open('data/bz%05d.d' %(it) ,'rb')
bz =np.fromfile(f, dtype='float64').reshape(-1,nx)

plt.imshow(bz,origin='lower',cmap='jet',aspect='auto',extent=([0,lx,0,ly]))
plt.xlabel(r'$x/(c/\omega_{pe})$'); plt.ylabel(r'$x/(c/\omega_{pe})$')
plt.title('$B_z$')
plt.colorbar()
plt.axis('square')
plt.show()

Snapshot in the Fourier space

The solid black line describes the theoretical curve of a numerical Cherenkov emission.
That is defined as follows: $${k_y}\Delta y = 2\arcsin \left\{\pm {\frac{{\Delta y}}{{c\Delta t}}\sqrt {{{\tan }^2}\left( {\frac{{V{k_x}\Delta t}}{2}} \right) - {{\left( {\frac{{c\Delta t}}{{\Delta x}}\sin \left( {\frac{{{k_x}\Delta x}}{2}} \right)} \right)}^2} - \frac{1}{2}\frac{{\omega _{pe}^2}}{{V{k_x}}}\Delta t\tan \left( {\frac{{V{k_x}\Delta t}}{2}} \right)} } \right\}$$

kmin=2*mt.pi/(dx*nx)*(-nx/2); kmax=2*mt.pi/(dx*nx)*(nx/2)
kx=np.linspace(kmin,kmax,nx)
a1=np.tan(kx*dt/2)**2
a2=-(dt/dx*np.sin(kx*dx/2))**2
a3=-0.5/kx*dt*np.tan(kx*dt/2)
kym=2*np.arcsin(+dy/dt*np.sqrt(a1+a2+a3))/dy
kyp=2*np.arcsin(-dy/dt*np.sqrt(a1+a2+a3))/dy

/home/mnakanot/anaconda3/lib/python3.6/site-packages/ipykernel_launcher.py:6: RuntimeWarning: invalid value encountered in sqrt
  
/home/mnakanot/anaconda3/lib/python3.6/site-packages/ipykernel_launcher.py:6: RuntimeWarning: invalid value encountered in arcsin
  
/home/mnakanot/anaconda3/lib/python3.6/site-packages/ipykernel_launcher.py:7: RuntimeWarning: invalid value encountered in sqrt
  import sys
/home/mnakanot/anaconda3/lib/python3.6/site-packages/ipykernel_launcher.py:7: RuntimeWarning: invalid value encountered in arcsin
  import sys

fftbz=fftshift(fft2(bz))
plt.imshow(abs(fftbz),aspect='auto',extent=([kmin,kmax,kmin,kmax]),cmap='jet',norm=mpl.colors.LogNorm())
plt.xlabel('$k_x(c/\omega_{pe})$');plt.ylabel('$k_y(c/\omega_{pe})$')
plt.colorbar()
plt.clim(1e-2,1e3)
plt.plot(kx,kym,'-k')
plt.plot(kx,kyp,'-k')
plt.show()

Time evolution of the magnetic field

tt=dt*np.arange(nd)
enbz=np.zeros((nd))
for it in range(nd):
    f =open('data/bz%05d.d' %(it) ,'rb')
    bz =np.fromfile(f, dtype='float64').reshape(-1,nx)
    enbz[it]=np.mean(bz**2)

plt.plot(tt,enbz)
plt.xlabel(r'$t\omega_e$'); plt.ylabel('$B_z$ energy')
plt.yscale('log')
plt.show()