M. Nakanotani
We consider electromagnetic ion cyclotron instabilities driven by a proton beam. The energy of the proton beam is transfered through waves which are resonant with the beam. The resonant wave grows in amplitudes and diffuses the proton beam. Here, we focus on 1. the linear growth of the ion cyclotron instability, 2. pitch angle scattering of a proton beam, and 3. nonlinear evolution of large-amplitude waves (decay instability).
We use 1D full PIC simulations with the following parameters.
The solid line is the dispersion relation of right- and left-handed circularly polarized waves, and dashed line represents the proton beam. The two intersections are a resonant condition of the proton beam and waves. At the intersection on the whistler branch, the proton beam energy is transfered to the wave mode (consequently, a wave excitation).
%config InlineBackend.figure_format='retina'
import math as mt
import matplotlib.pyplot as plt
import numpy as np
import scipy.fft as sf
from scipy.fftpack import fft, ifft,fft2,fftshift,ifft2,ifftshift
import datetime
import bisect
import matplotlib.patches as patches
import matplotlib.dates as mdates
from IPython.display import display,Math
import matplotlib as mpl
import pandas as pd
from matplotlib.gridspec import GridSpec
import matplotlib.animation as animation
from mpl_toolkits.axes_grid1.inset_locator import zoomed_inset_axes, mark_inset
from matplotlib.ticker import ScalarFormatter
from matplotlib.colors import ListedColormap, LinearSegmentedColormap
import plasmapy.dispersion.dispersionfunction as pdd
from scipy.special import spherical_in
import scipy
from scipy.integrate import quad
plt.rcParams['font.size'] = 14
plt.rcParams['axes.linewidth'] = 1.5
plt.rcParams['animation.embed_limit']=60
plt.rcParams['animation.html'] = 'jshtml'
kk=np.linspace(-5,5,1000)
wwp=np.zeros(1000)
wwm=np.zeros(1000)
for ik in range(1000):
k2=kk[ik]**2
wwp[ik]=(k2+np.sqrt(k2**2+4*k2))/2
wwm[ik]=(k2-np.sqrt(k2**2+4*k2))/2
fig = plt.figure()
ax = fig.add_subplot(1, 1, 1)
ax.plot(kk,wwp,'k')
ax.plot(kk,wwm,'k')
ax.plot(kk,3*kk-1,'k--')
#ax.set_xlim(0,2)
ax.set_ylim(-2,4)
ax.spines['left'].set_position('center')
ax.spines['bottom'].set_position('zero')
ax.spines['right'].set_color('none')
ax.spines['top'].set_color('none')
#plt.tick_params(left = False, right = False , labelleft = False ,
# labelbottom = False, bottom = False)#ax.yaxis.set_ticks_position('left')
ax.text(5.5,-0.4,'k')
ax.text(0.1,4.2,'$\omega$')
plt.show()
The following plot is the dispersion relation of the unstable mode explained above. This can be written as follows,
$$ \omega^2-k^2c^2+\frac{n_b}{n_0}\omega_{pi}^2\frac{\omega-kV_b}{kV_{thb}}Z\left(\frac{\omega-kV_b+\Omega_i}{kV_{thb}}\right)+\left(1-\frac{n_b}{n_0}\right)\omega^2_{pi}\frac{\omega}{kV_{thi}}Z\left(\frac{\omega+\Omega_i}{kV_{thi}}\right)+\omega^2_{pe}\frac{\omega-kV_{ed}}{kV_{the}}Z\left(\frac{\omega-kV_{ed}-\Omega_e}{kV_{the}}\right)=0. $$As you can see, the maximum growth rate is $0.4\Omega_{i}$ at $k=0.8(v_A/\Omega_i)$.
def disp(W,k):
c=1
xi=(W-k*vi)/(k*vti)
Zi=pdd.plasma_dispersion_func(xi+wci/(k*vti))
Di=(wpi/W)**2*xi*Zi
xb=(W-k*vb)/(k*vtb)
Zb=pdd.plasma_dispersion_func(xb+wcb/(k*vtb))
Db=(wpb/W)**2*xb*Zb
xe=(W-k*ve)/(k*vte)
Ze=pdd.plasma_dispersion_func(xe+wce/(k*vte))
De=(wpe/W)**2*xe*Ze
return 1-(k*c/W)**2+(Di+Db+De)
from scipy import optimize
b0=0.2
rm=25
nb=0.2
bte=0.5
bti=0.03125
vA =b0/np.sqrt(rm)
vte=np.sqrt(2)*np.sqrt(bte/2)*b0
vti=np.sqrt(2)*np.sqrt(bti/2)*b0/np.sqrt(rm)
vtb=vti
wpe=1
wpi= wpe/np.sqrt(rm)*np.sqrt(1-nb)
wpb= wpe/np.sqrt(rm)*np.sqrt(nb)
wce=-wpe*b0
wci= wpe*b0/rm
wcb= wpe*b0/rm
vb=3*vA
vi=0#3*vA*nb
ve=vb*nb
k0=0.6/(vA/wci)
w0=0.6*wci
kk1=np.linspace(k0,4/(vA/wci),200)
ww1=np.zeros(200,dtype=complex)
kk2=np.linspace(k0,2e-2/(vA/wci),100)
ww2=np.zeros(100,dtype=complex)
sol=optimize.newton(lambda x: disp(x,k0), x0=w0)
print(sol)
for ik in range(len(kk1)):
sol=optimize.newton(lambda x: disp(x,kk1[ik]), x0=sol, maxiter=100)
ww1[ik]=sol
sol=optimize.newton(lambda x: disp(x,k0), x0=w0)
for ik in range(len(kk2)):
sol=optimize.newton(lambda x: disp(x,kk2[ik]), x0=sol,maxiter=100)
ww2[ik]=sol
kk=np.concatenate([kk2[::-1],kk1])
ww=np.concatenate([ww2[::-1],ww1])
idx=np.argmax(((np.imag(ww))))
plt.plot(kk*(vA/wci),np.real(ww)/wci,label='$\omega_r$')
plt.plot(kk*(vA/wci),np.imag(ww)/wci,label='$\gamma$')
#plt.plot(kk*(vA/wci),-1+3*kk*(vA/wci))
plt.legend()
plt.xlabel(r'$k(v_A/\Omega_i)$')
plt.ylabel(r'$\omega/\Omega_i$')
plt.xlim(0,3)
plt.ylim(-0,0.5)
plt.grid()
plt.show()
nx=1024
nt=80000
isavw=100
isavp=1000
dx=0.1
dt=0.1
rm=25
b0=0.2
def anim_phase(it):
f1=open('data/up1%05d.d' %(it),'rb')
f2=open('data/up2%05d.d' %(it),'rb')
f3=open('data/up3%05d.d' %(it),'rb')
xp1=np.fromfile(f1,dtype='float64').reshape(-1,5)
xp2=np.fromfile(f2,dtype='float64').reshape(-1,5)
xp3=np.fromfile(f3,dtype='float64').reshape(-1,5)
fig.clf()
ax1 = fig.add_subplot(111)
ax1.clear()
xx=dx*np.arange(nx)
ax1.plot(xp2[:,0],xp2[:,1],',',label='core proton')
ax1.plot(xp3[:,0],xp3[:,1],',',label='beam proton')
ax1.legend()
ax1.plot(xx,byhst[it*10,:],'-')
ax1.plot(xx,bzhst[it*10,:],'-')
ax1.set_xlim(0,nx*dx)
ax1.set_ylim(-0.25,0.25)
ax1.set_xlabel('$X$')
ax1.set_ylabel('$V_x$')
f5 =open('data/byhst.d' ,'rb')
f6 =open('data/bzhst.d' ,'rb')
byhst =np.fromfile(f5, dtype='float64' ).reshape(-1,nx)
bzhst =np.fromfile(f6, dtype='float64' ).reshape(-1,nx)
fig=plt.figure(figsize=(16,6))
anim=animation.FuncAnimation(fig,anim_phase,frames=80)
plt.close()
anim
The energy plot below for each component shows that the proton beam energy is primary transfered to magnetic fields and core protons. The linear stage of the instability is in the range of $2000<time<2500$. After the linear stage, the electron energy gradually increases. The relative error stays within $0.3$%.
tt=(dt*isavw)*np.arange(nt//isavw+1)
f1=open('data/enhst.d','rb')
enhst=np.fromfile(f1,dtype='float64').reshape(-1,5)
plt.figure(figsize=(12,6))
plt.plot(tt,enhst[:,0],label='electron')
plt.plot(tt,enhst[:,1],label='core proton')
plt.plot(tt,enhst[:,2],label='beam proton')
plt.plot(tt,enhst[:,3],label='electric field')
plt.plot(tt,enhst[:,4],label='magnetic field')
plt.xlabel('time')
plt.ylabel('energy')
plt.legend()
plt.show()
entot=np.sum(enhst,axis=1)
plt.figure(figsize=(12,6))
plt.plot(tt,(1-entot/entot[1])*100,'-')
plt.xlabel('time')
plt.ylabel('relative error [%]')
plt.show()
Linear growth calculated from the simulation agrees well with the theoretical linear growth.
f1 =open('data/exhst.d' ,'rb')
f2 =open('data/eyhst.d' ,'rb')
f3 =open('data/ezhst.d' ,'rb')
f4 =open('data/bxhst.d' ,'rb')
f5 =open('data/byhst.d' ,'rb')
f6 =open('data/bzhst.d' ,'rb')
f7 =open('data/rohst.d' ,'rb')
f8 =open('data/jxhst.d' ,'rb')
f9 =open('data/jyhst.d' ,'rb')
f10=open('data/jzhst.d' ,'rb')
exhst =np.fromfile(f1, dtype='float64' ).reshape(-1,nx)
eyhst =np.fromfile(f2, dtype='float64' ).reshape(-1,nx)
ezhst =np.fromfile(f3, dtype='float64' ).reshape(-1,nx)
bxhst =np.fromfile(f4, dtype='float64' ).reshape(-1,nx)
byhst =np.fromfile(f5, dtype='float64' ).reshape(-1,nx)
bzhst =np.fromfile(f6, dtype='float64' ).reshape(-1,nx)
rohst =np.fromfile(f7, dtype='float64' ).reshape(-1,nx)
jxhst =np.fromfile(f8, dtype='float64' ).reshape(-1,nx)
jyhst =np.fromfile(f9, dtype='float64' ).reshape(-1,nx)
jzhst =np.fromfile(f10,dtype='float64' ).reshape(-1,nx)
### Fourier decomposition (Terasawa1986JGR)
import scipy.integrate as integrate
import scipy.special as special
lx=nx*dx
x=np.linspace(-lx/2,lx/2,nx)
byc=np.zeros((nt//isavw,nx//2))
bys=np.zeros((nt//isavw,nx//2))
bzc=np.zeros((nt//isavw,nx//2))
bzs=np.zeros((nt//isavw,nx//2))
for ik in range(nx//2):
byc[:,ik]=2/lx*np.sum(byhst[0:nt//isavw,:]*np.cos(2*np.pi*(ik+1)*x/lx)*dx,axis=1)
bys[:,ik]=2/lx*np.sum(byhst[0:nt//isavw,:]*np.sin(2*np.pi*(ik+1)*x/lx)*dx,axis=1)
bzc[:,ik]=2/lx*np.sum(bzhst[0:nt//isavw,:]*np.cos(2*np.pi*(ik+1)*x/lx)*dx,axis=1)
bzs[:,ik]=2/lx*np.sum(bzhst[0:nt//isavw,:]*np.sin(2*np.pi*(ik+1)*x/lx)*dx,axis=1)
fftbr=0.5*((byc+bzs)+1j*(bzc-bys))
fftbl=0.5*((byc-bzs)+1j*(bzc+bys))
br=np.zeros((nt//isavw,nx),dtype='complex')
bl=np.zeros((nt//isavw,nx),dtype='complex')
for it in range(nt//isavw):
for ik in range(nx//2):
br[it,:]=br[it,:]+fftbr[it,ik]*np.exp( 1j*2*np.pi*(ik+1)*x/lx)
bl[it,:]=bl[it,:]+fftbl[it,ik]*np.exp(-1j*2*np.pi*(ik+1)*x/lx)
byr=np.real(br)
bzr=np.imag(br)
byl=np.real(bl)
bzl=np.imag(bl)
wci=b0/rm
kmin=2*mt.pi/(dx*nx)*(-nx/2);kmax=2*mt.pi/(dx*nx)*(nx/2)
kmin=kmin*np.sqrt(rm); kmax=kmax*np.sqrt(rm)
kk=np.linspace(kmin,kmax,nx)
tt=(dt*isavw)*np.arange(nt//isavw)
plt.figure(figsize=(5,3))
#plt.rcParams['font.size'] = 16
fftbzsave=fftshift(fft(byl),axes=1)
#idx=np.argmax(np.imag(roots[:,6]))#
idx=np.abs(kk-1).argmin()
plt.figure(figsize=(12,6))
#for i in range(20):
# plt.semilogy(tt,abs(fftbzsave[:,2050+i])*(4**i),'-')
plt.semilogy(tt,abs(fftbzsave[:,idx]),'-')
#plt.semilogy(tt,np.mean(byl**2,axis=1))
plt.xlabel(r'$time$')
plt.ylabel(r'$|B_y^{FFT}|$')
plt.plot(tt,2e-2*np.exp(0.4*wci*tt),'--k',lw=1,label=r'$\gamma\sim0.4\Omega_i$')
plt.ylim(1e-2,1e3)
#plt.axis([0,20*dt,1e-2,3e0])
plt.legend(frameon=False)
#plt.xlim(0,4)
#plt.title("Time evolution of the amplitude of the maximum mode (solid line)")
#plt.savefig('figure/Ex_t.pdf',format='pdf',dpi=200,bbox_inches='tight')
plt.show()
The hodograph of $B_y-B_z$ verifies that the excited wave is circularly right-handed polarized.
plt.plot(byhst[0:300,0],bzhst[0:300,0])
plt.xlim(-0.2,0.2)
plt.ylim(-0.2,0.2)
plt.xlabel('$B_y$')
plt.ylabel('$B_z$')
plt.show()
Due to the excited waves, the proton beam is scattered in the phase space.
def anim_pitch(it):
f1=open('data/up1%05d.d' %(it),'rb')
f2=open('data/up2%05d.d' %(it),'rb')
f3=open('data/up3%05d.d' %(it),'rb')
xp1=np.fromfile(f1,dtype='float64').reshape(-1,5)
xp2=np.fromfile(f2,dtype='float64').reshape(-1,5)
xp3=np.fromfile(f3,dtype='float64').reshape(-1,5)
fig.clf()
ax1 = fig.add_subplot(111)
ax1.clear()
xion=np.r_[xp2,xp3]
weight=np.r_[0.8*np.ones(len(xp2)),0.2*np.ones(len(xp3))]
ax1.hist2d(xion[:,1],xion[:,2],bins=256,weights=weight,cmap='terrain_r',norm=mpl.colors.LogNorm())
th=np.linspace(0,2*np.pi,100)
for ir in range(18):
ax1.plot(ir*0.01*(np.cos(th))+0.07 ,ir*0.01*np.sin(th),'k',lw=0.5)
ax1.set_xlabel(r'$v_x/c$')
ax1.set_ylabel(r'$v_y/c$')
ax1.set_xlim(-0.08,0.18)
ax1.set_ylim(-0.08,0.08)
fig=plt.figure(figsize=(10,5))
anim=animation.FuncAnimation(fig,anim_pitch,frames=120)
plt.close()
anim