I reproduce a subject in this page http://www.astro.phys.s.chiba-u.ac.jp/pcans/application_em1d.html.
I use python3.
%config InlineBackend.figure_format = 'retina'
from numba import jit,f8,u8
import numpy as Np
import math as mt
import matplotlib.pyplot as plt
import math as mt
from scipy.fftpack import fft, ifft,fft2,fftshift
import time
from tqdm import tqdm
import matplotlib.animation as animation
#Electromagnetic 1-d PIC code
#@jit#((f8,f8,f8,f8,u8,u8,u8))
def em1(rm,b0,th,vta,vte,ppc,nx,nt):
Np.seterr(divide='ignore', invalid='ignore')
#-----
#---wpe=1, c=1, me=1
#-----
np=ppc*nx
dx=vta; dt=dx;lx=dx*nx
bx0=b0*mt.cos(th*mt.pi/180); bz0=b0*mt.sin(th*mt.pi/180)
xx=dx*Np.arange(nx); tt=dt*Np.arange(nt)
xp,vx,vy,vz=Np.zeros((4,2,np))
xp[0,:]=Np.linspace(0,lx-lx/np,np)
vx[0,:]=Np.random.normal(0,vta,np)
vy[0,:]=Np.random.normal(0,vte,np)
vz[0,:]=Np.random.normal(0,vte,np)
xp[1,:]=Np.linspace(0,lx-lx/np,np)
vx[1,:]=Np.random.normal(0,vta,np)/mt.sqrt(rm)
vy[1,:]=Np.random.normal(0,vta,np)/mt.sqrt(rm)
vz[1,:]=Np.random.normal(0,vta,np)/mt.sqrt(rm)
gm=Np.sqrt(1+vx**2+vy**2+vz**2)
vx=vx/gm; vy=vy/gm; vz=vz/gm
jj=Np.arange(nx); jp=Np.r_[Np.arange(nx-1)+1, 0]; jm=Np.r_[nx-1, Np.arange(nx-1)]
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
q,qomdt,qdt=Np.zeros((3,2))
q[0]=-nx/np; qomdt[0]=-0.5*dt; qdt[0]=0.25*dt*q[0]
q[1]=-q[0]; qomdt[1]= 0.5*dt/rm; qdt[1]=0.25*dt*q[1]
ex,ey,ez=Np.zeros((3,nx))
by,bz =Np.zeros((2,nx))
eym,eyp =Np.zeros((2,nx))
ezm,ezp =Np.zeros((2,nx))
exsave,eysave,ezsave =Np.zeros((3,nt,nx))
bzsave,bysave,enesave=Np.zeros((3,nt,nx))
tasave,tesave=Np.zeros((2,nt))
dssave=Np.zeros((nt,nx))
print('initial condition')
plt.subplot(121);plt.plot(xp[0,:],vx[0,:],',');plt.plot(xp[0,:],vy[0,:],',')
plt.subplot(122);plt.hist(vx[0,:],bins=500,histtype='step');plt.hist(vy[0,:],bins=500,histtype='step')
plt.show()
for it in tqdm(range(nt)):
#---solve eqs. of motion for each particle
vx,vy,vz,gm=acc(np,nx,dx,qomdt,bx0,bz0,xp,vx,vy,vz,gm,ex,ey,ez,by,bz)
#enesave[it,5+isp]=0.5*enesave[it,5+isp]/np
#---calculate current on each grid
jym,jzm=Np.zeros((2,2,nx))
jym,jzm=currnt(np,nx,dx,qdt,xp,vy,vz,jym,jzm)
#---calculate new particle posotions
xp=xp+dt*vx
xp[xp>lx]=xp[xp>lx]-lx
xp[xp<0 ]=xp[xp<0 ]+lx
##---calculate charge density and current---#
ds =Np.zeros((2,nx))
jyp,jzp=Np.zeros((2,2,nx))
ds,jyp,jzp=currntden(np,nx,dx,q,qdt,xp,vy,vz,jyp,jzp,ds)
#---solve Maxwell's equations on each grid
#---longitudinal
#exfft=-1j/kk*fft(Np.sum(ds,0))
exfft=-1j*kl*fft(Np.sum(ds,0))/Kl
exfft[0]=0
ex=Np.real(ifft(exfft))
exsave[it,:]=ex
#---transverse
eym=eym-Np.sum(jym,0)
eyp=eyp-Np.sum(jym,0)
ezm=ezm-Np.sum(jzm,0)
ezp=ezp-Np.sum(jzm,0)
eym[jj]=eym[jp]
ezm[jj]=ezm[jp]
eyp[jj]=eyp[jm]
ezp[jj]=ezp[jm]
eym=eym-Np.sum(jyp,0)
eyp=eyp-Np.sum(jyp,0)
ezm=ezm-Np.sum(jzp,0)
ezp=ezp-Np.sum(jzp,0)
ey=eyp+eym
bz=eyp-eym
ez=ezp+ezm
by=ezm-ezp
eysave[it,:]=ey
ezsave[it,:]=ez
bysave[it,:]=by
bzsave[it,:]=bz
#---calculate electron temperature parallel and perpendicular and density
tasave[it]=Np.average(vx[0,:]**2)
tesave[it]=Np.average((vy[0,:]**2+vz[0,:]**2)/2)
dssave[it,:]=ds[0,:]
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)
kaxis=Np.linspace(kmin,kmax,nx); waxis=Np.linspace(wmin,wmax,nt)
return locals()
@jit('f8[:,:](u8,u8,f8,f8[:],f8,f8,f8[:,:],f8[:,:],f8[:,:],f8[:,:],f8[:,:],f8[:],f8[:],f8[:],f8[:],f8[:])')
def acc(np,nx,dx,qomdt,bx0,bz0,xp,vx,vy,vz,gm,ex,ey,ez,by,bz):
for isp in range(2):
for ip in range(np):
ixm=int(xp[isp,ip]/dx); ixp=ixm+1
wxp=xp[isp,ip]/dx-ixm; wxm=1-wxp
if ixp>nx-1: ixp=ixp-nx
exap=qomdt[isp]*(wxm*ex[ixm]+wxp*ex[ixp])
eyap=qomdt[isp]*(wxm*ey[ixm]+wxp*ey[ixp])
ezap=qomdt[isp]*(wxm*ez[ixm]+wxp*ez[ixp])
bxap=qomdt[isp]*bx0
byap=qomdt[isp]*(wxm*by[ixm]+wxp*by[ixp])
bzap=qomdt[isp]*(wxm*bz[ixm]+wxp*bz[ixp]+bz0)
#---half acceleration
gvxs=vx[isp,ip]*gm[isp,ip]+exap
gvys=vy[isp,ip]*gm[isp,ip]+eyap
gvzs=vz[isp,ip]*gm[isp,ip]+ezap
#---first half rotation
gmm=mt.sqrt(1+gvxs**2+gvys**2+gvzs**2)
#enesave[it,5+isp]=enesave[it,5+isp]+gm[isp,ip]-1
bxap=bxap/gmm
byap=byap/gmm
bzap=bzap/gmm
vvx=gvxs+gvys*bzap-gvzs*byap
vvy=gvys+gvzs*bxap-gvxs*bzap
vvz=gvzs+gvxs*byap-gvys*bxap
#---second half rotation + second half acceleration
f=2/(1+bxap**2+byap**2+bzap**2)
bxap=f*bxap
byap=f*byap
bzap=f*bzap
gvxs=gvxs+vvy*bzap-vvz*byap+exap
gvys=gvys+vvz*bxap-vvx*bzap+eyap
gvzs=gvzs+vvx*byap-vvy*bxap+ezap
#---half acceleration
gm[isp,ip]=mt.sqrt(1+gvxs**2+gvys**2+gvzs**2)
vx[isp,ip]=gvxs/gm[isp,ip]
vy[isp,ip]=gvys/gm[isp,ip]
vz[isp,ip]=gvzs/gm[isp,ip]
return vx,vy,vz,gm
@jit('f8[:,:](u8,u8,f8,f8[:],f8[:,:],f8[:,:],f8[:,:],f8[:,:],f8[:,:])')
def currnt(np,nx,dx,qdt,xp,vy,vz,jym,jzm):
for isp in range(2):
for ip in range(np):
ixm=mt.floor(xp[isp,ip]/dx); ixp=ixm+1
#print(ixm,ixp)
wxp=xp[isp,ip]/dx-ixm; wxm=1-wxp
if ixp>nx-1: ixp=ixp-nx#; print(ixp)
jym[isp,ixm]=jym[isp,ixm]+wxm*vy[isp,ip]
jym[isp,ixp]=jym[isp,ixp]+wxp*vy[isp,ip]
jzm[isp,ixm]=jzm[isp,ixm]+wxm*vz[isp,ip]
jzm[isp,ixp]=jzm[isp,ixp]+wxp*vz[isp,ip]
jym[isp,:]=jym[isp,:]*qdt[isp]
jzm[isp,:]=jzm[isp,:]*qdt[isp]
return jym,jzm
@jit('f8[:,:](u8,u8,f8,f8[:],f8[:],f8[:,:],f8[:,:],f8[:,:],f8[:,:],f8[:,:],f8[:,:])')
def currntden(np,nx,dx,q,qdt,xp,vy,vz,jyp,jzp,ds):
for isp in range(2):
for ip in range(np):
#print(xp[ip])
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#; print(ixp)
ds[isp,ixm]=ds[isp,ixm]+wxm
ds[isp,ixp]=ds[isp,ixp]+wxp
jyp[isp,ixm]=jyp[isp,ixm]+wxm*vy[isp,ip]
jyp[isp,ixp]=jyp[isp,ixp]+wxp*vy[isp,ip]
jzp[isp,ixm]=jzp[isp,ixm]+wxm*vz[isp,ip]
jzp[isp,ixp]=jzp[isp,ixp]+wxp*vz[isp,ip]
ds[isp,:] =ds[isp,:]*q[isp]
jyp[isp,:]=jyp[isp,:]*qdt[isp]
jzp[isp,:]=jzp[isp,:]*qdt[isp]
return ds,jyp,jzp
Yoshizumi Miyoshi(Nagoya University) and Takashi Minoshima(JAMSTEC)
http://www.astro.phys.s.chiba-u.ac.jp/pcans/em1d_whistler.html
$$f(v_{\perp},v_{||})=2\pi n\left(\sqrt{\frac{m}{2\pi\kappa T_\perp}}\right)^2\sqrt{\frac{m}{2\pi\kappa T_{||}}} v_\perp\exp\left(-\frac{mv_\perp^2}{2\kappa T_\perp}-\frac{mv_{||}^2}{2\kappa T_{||}}\right)$$
$$\rm{First:}\quad\frac{T_{\perp}}{B}\sim\rm{constant}$$
$$\rm{Second:}\quad\frac{T_{||}B^2}{n^2}\sim\rm{constant}$$
The linear growth rate of the whistler waves with an electron temperature anisotropy can be obtained by solving Vlasov equation and Maxwell equations (Kennel & Petschek, 1966).
$$\omega_r=\frac{k^2c^2}{\omega_e^2+k^2c^2}\Omega_e, \quad \omega_i=\pi\Omega_e\left(1-\frac{\omega_r}{\Omega_e}\right)^2\eta(V_R)[A(V_R)-A_c].$$
Here $\omega_e$ is the electron plasma frequency and $\Omega_e$ is the electron cyclotron frequency.
$$\eta ({V_R}) = {\left. {2\pi \frac{{{\Omega _e} - {\omega _r}}}{k}\int_0^\infty {{v_ \bot }} Fd{v_ \bot }} \right|_{{v_{||}} = {V_R}}}, \quad A({V_R}) = {\left. {\frac{{\int_0^\infty {\left( {{v_{||}}\frac{{\partial F}}{{\partial {v_ \bot }}} - {v_ \bot }\frac{{\partial F}}{{\partial {v_{||}}}}} \right)\frac{{v_ \bot ^2}}{{{v_{||}}}}d{v_ \bot }} }}{{2\int_0^\infty {{v_ \bot }Fd{v_ \bot }} }}} \right|_{{v_{||}} = {V_R}}},\quad A_c=\frac{1}{\Omega_e/\omega_r-1}.$$
The resonance speed $V_R$ satisfies $kV_R=\omega_r-\Omega_e$.
If the distribution function $F$ is Maxwell distribution, $A$ can be written as $A=T_\perp/T_{||}-1$.
The whistler mode is in a root of the dispersion relation described below (${\bf k}=(0,0,k), {\bf B}_0=(0,0,B_0)$): $$\left(\frac{ck}{\omega}\right)^2=R=1-\sum_s\frac{\omega_{ps}^2}{\omega(\omega+\Omega_s)},\quad {\bf\tilde{E}}=(E_0,{\rm i}E_0,0).$$ We can approximate the equation if $$\Omega_i\ll\omega<\Omega_e\quad \rm{and} \quad \frac{\omega_{pe}}{\Omega_e}\gg1.$$ It becomes $$\omega=\frac{k^2c^2}{\omega_{pe}^2+k^2c^2}\Omega_e.$$
We show an example of the dispersion relation of the linear growth rate.
Here the ration of the plasma frequency and the cyclotron frequency is 5 and the temperature anisotropy $A=1.2$.
import matplotlib.pyplot as plt
import numpy as np
import math as mt
vta=0.1*mt.sqrt(2)
A=1.2
k =np.linspace(0.01,1.75,1000)
wr=k**2/(1+k**2)
VR=(1-wr)/k/5
eta=1/np.sqrt(np.pi)*VR/vta*np.exp(-VR**2/vta**2)
wi=mt.pi*(1-wr)**2*eta*(A-1/(1/wr-1))
plt.figure(figsize=(12,5))
plt.rcParams['font.size'] = 16
plt.subplot(211); plt.plot(k,wr,'-k')
plt.ylabel(r'$\omega_r/\Omega_e$')
plt.ylim(-0.,1)
plt.grid()
plt.subplot(212); plt.plot(k,wi,'-k')
plt.xlabel(r'$k(c/\omega_{pe})$')
plt.ylabel(r'$\omega_i/\Omega_e$')
plt.ylim(-0.,1)
plt.grid()
plt.show()
As the initial condition, we use the paramters below referring Sydora et al. (2007).
rm=1840; b0=0.2; th=0.0; vta=0.14; vte=0.21; ppc=300; nx=1024; nt=4096
%time data=em1(rm,b0,th,vta,vte,ppc,nx,nt) #rm,b0,th,vta,vte,ppc,nx,nt
locals().update(data)
#---make an animation---#
fig=plt.figure(figsize=(15,5))
ax = fig.add_subplot(111)
def update_anim(i):
it=16*i
ax.cla()
ax.plot(xx,bysave[it,:],'b-')
ax.plot(xx,bzsave[it,:],'r-')
ax.plot(xx, Np.sqrt(bysave[it,:]**2+bzsave[it,:]**2),'k-')
ax.plot(xx,-Np.sqrt(bysave[it,:]**2+bzsave[it,:]**2),'k-')
ax.set_ylim(-0.1,0.1)
ax.set_xlabel(r'$x/(c/\omega_{pe})$')
ax.set_ylabel(r'$B_y, B_z$')
plt.rcParams['animation.html'] = 'html5'
anim = animation.FuncAnimation(fig, update_anim, blit=False,frames=256)
plt.close()
#ims=[]
#fig=plt.figure(figsize=(15,5))
#for it in range(nt):
# if(it%16==0):
# im=plt.plot(xx,bysave[it,:],'k-')
# plt.xlabel(r'$x/(c/\omega_{pe})$')
# plt.ylabel(r'$B_y$')
# ims.append(im)
#
#plt.rcParams['animation.html'] = 'html5'
#ani=animation.ArtistAnimation(fig,ims)
#plt.close()
anim
plt.figure(figsize=(12,10))
plt.imshow(bzsave,interpolation='none',extent=[xx[0],xx[nx-1],tt[0],tt[nt-1]/5],
origin='lower',aspect='auto',cmap='jet')
plt.axis([0, xx[nx-1], 0, tt[nt-1]/5])
plt.xlabel(r'$x/(c/\omega_{pe})$')
plt.ylabel(r'$T\Omega_{e}$')
plt.title('Time-space evolution of $B_y$')
plt.colorbar()
plt.show()
import numpy as np
plt.figure(figsize=(12,5))
plt.rcParams['font.size'] = 16
fftbzsave=fftshift(fft(bysave),axes=(1,))
idx=np.abs(kk-0.74).argmin()
plt.semilogy(tt/5,abs(fftbzsave[:,512+idx])/0.2,'-')
plt.xlabel(r'$T\Omega_{e}$')
plt.ylabel(r'$B_{yk}/B_0$')
plt.plot(tt/5,np.exp(tt/5*0.18),'--k')
plt.axis([0,60,0,100])
plt.title("Time evolution of the amplitude of the maximum mode (solid line)")
plt.show()
From the panel, we can see that the simulation result and the theory is well consisted at the linear phase.
plt.figure(figsize=(12,6))
plt.rcParams['font.size']=16
plt.title(r"The time evolution of $T_\perp/T_{||}$")
plt.plot(tt/5,tesave/tasave,'-k')
plt.ylabel(r'$T_{\perp}/T_{||}$')
plt.title(r"The time evolution of $T_\perp/T_{||}$")
#plt.show()
#plt.plot(tt/5,tesave,'-k',label=r'$T_\perp$')
#plt.plot(tt/5,tasave,'-b',label=r'$T_{||}$')
#plt.legend()
#plt.ylabel(r'$T_{\perp},T_{||}$')
plt.xlabel(r'$T\Omega_{e}$')
plt.show()
Free energy to drive waves is the temperature anisotropy of a distribution function.
From the panel, as the waves grow, the ration of the perpendicular and parallel temeprature is being relaxed.
The ration is about 1.2 at the saturation stage.
fftbzsave=fftshift(fft(bzsave),axes=(1,))
#plt.plot(tt,abs(fftbzsave[:,256]))
#plt.pcolor(abs(fftbzsave))
plt.figure(figsize=(12,6))
plt.rcParams['font.size']=16
plt.imshow(np.log(0.01+abs(fftbzsave)),interpolation='none',extent=[kmin,kmax,tt[0],tt[nt-1]/5],
origin='lower',aspect='auto',cmap='jet')
plt.axis([0, kmax/15, 0, tt[nt-1]/5])
plt.colorbar()
plt.clim(-4,3)
plt.xlabel(r'$k(c/\omega_{pe})$')
plt.ylabel(r'$T\Omega_{e}$')
plt.title("The time evolution of the wave spectral")
plt.show()
In the linear phase, waves are drived across wavenumber band focusing on $kc/\omega_{pe}\sim 0.7$.
In the nonlinear phase, they shift to long-wavelength mode (inverse cascade).
The reason why is that as the anisotropy is relaxed, the maximum growth wavelength shifts to long-wavelength mode (Sydora et al., 2007).
fftbzsave=fftshift(fft2(bzsave))
plt.figure(figsize=(12,6))
plt.rcParams['font.size']=16
plt.imshow(np.log(0.01+abs(fftbzsave)),interpolation='none',extent=[kmin,kmax,wmin*5,wmax*5],
aspect='auto',cmap='jet')
plt.axis([0, kmax/15, 0, wmax/20])
k =np.linspace(0.01,1.75,1000)
wr=k**2/(1+k**2)
plt.plot(k,wr,'-k')
plt.plot(k,wi,'--k')
plt.xlabel(r'$k(c/\omega_{pe})$')
plt.ylabel(r'$\omega/\Omega_{e}$')
plt.clim(1,10)
#plt.colorbar()
plt.title("A dispersion relation of waves driven by the electron temerpature anisotropy")
plt.show()
The solid line shows R-mode (whistler wave) which is one of low-frequency electromagnetic modes.
We can see that waves driven by electron temperature anisotropy is whislter waves.
http://slideplayer.com/slide/10107285/
http://www.igpp.ucla.edu/public/mkivelso/refs/PUBLICATIONS/kennelPetschekJZ071i001p00001.pdf
http://www.johnstonsarchive.net/physics/RBtutorial.pdf
http://slideplayer.com/slide/7543046/
http://rdecw.nifs.ac.jp/frontier/2012/2012files/miyoshi2012.pdf
http://www.spacestorm.eu/wp-content/uploads/sites/26/2016/06/Pres_46_Horne_IOP_20151.pdf
https://www.ncbi.nlm.nih.gov/pmc/articles/PMC4497482/pdf/jgra0119-1587.pdf
http://onlinelibrary.wiley.com/book/10.1029/GM169
http://onlinelibrary.wiley.com/doi/10.1029/2009JA014917/epdf
https://arxiv.org/pdf/1512.00919.pdf