I reproduce a subject in this page http://www.astro.phys.s.chiba-u.ac.jp/pcans/application_em1d.html.
I use python3.

In [1]:
%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
In [2]:
#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

Electron Temperature Anisotropy Instability

Yoshizumi Miyoshi(Nagoya University) and Takashi Minoshima(JAMSTEC)

http://www.astro.phys.s.chiba-u.ac.jp/pcans/em1d_whistler.html

Background

  • An temperature isotropy($T_{\perp}/T_{||}>1$) is oftenly observed in the Earth magnetoshpere.
  • Once charged particles are convected to stronger magnetic field, during keeping the first and second adiabatic invariances, the anisotropy is naturaly produced.
  • The anisotropy can drive ion cyclotron waves for ions and whistler waves for electrons.
  • Recently, frequency drifiting called "chorus" which frequencies of whistler waves drift in the nonlinear phase of the instability is atracts attention because of an acceleration of relativistic electrons by the chorus waves at the Earth radiation belt (Omura et al., 2008).

Bi-Maxwellian

$$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)$$

The First and Second Adiabatic Invariants

$$\rm{First:}\quad\frac{T_{\perp}}{B}\sim\rm{constant}$$

$$\rm{Second:}\quad\frac{T_{||}B^2}{n^2}\sim\rm{constant}$$

Linear Growth Rate

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$.

Waves Propagation Parallel to the Magnetic Field

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$.

In [3]:
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()

1D full-PIC simulation

As the initial condition, we use the paramters below referring Sydora et al. (2007).

  • $\Omega_{ce}/\omega_{pe}=0.2$
  • $T_\perp/T_{||}=2.2$
  • $T_{e||}/T_i=1$
  • $m_i/m_e=1840$
  • $\beta_i=\beta_{e||}=1.0$
In [4]:
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)

initial condition

100%|██████████| 4096/4096 [04:51<00:00, 14.05it/s]
CPU times: user 4min 23s, sys: 3.06 s, total: 4min 26s
Wall time: 4min 52s



















In [5]:



    


#---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()

















In [6]:



    


anim




















Out[6]:

























In [7]:



    


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()














































In [8]:



    


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.













In [9]:



    


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.













In [10]:



    


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).













In [11]:



    


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.

















References


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