Collisionless shock waves

I show a example of a collisionless shock simulation by using a 1d PIC simulation code.
The parameters are :

  • $\Omega_e/\omega_e=0.5$
  • $\Theta_{Bn}=90^\circ$
  • $m_i/m_e=400$
  • $v_{bulk}/c=0.1$
  • $\beta_i=\beta_e=0.02$
  • $v_A=0.025$
  • $M_A=5\sim6$
  • $ppc=20$.
In [1]:
%config InlineBackend.figure_format='retina'
import numpy as np
import matplotlib.pyplot as plt
from IPython import display
import matplotlib.animation as animation
Overview

We show the magnetic field component $B_z$ in space-time.
The shock reformation can be seen clearly.

In [2]:
%%time
nt=485
exsave=np.zeros((nt,20000))
eysave=np.zeros((nt,20000))
ezsave=np.zeros((nt,20000))
bysave=np.zeros((nt,20000))
bzsave=np.zeros((nt,20000))
for it in range(nt):
    f1=open('ex_%05d' %(it),'rb');f2=open('ey_%05d' %(it),'rb');f3=open('ez_%05d' %(it),'rb')
    f4=open('by_%05d' %(it),'rb');f5=open('bz_%05d' %(it),'rb')
    ex=np.fromfile(f1,dtype='float64').reshape(-1,2)
    ey=np.fromfile(f2,dtype='float64').reshape(-1,2)
    ez=np.fromfile(f3,dtype='float64').reshape(-1,2)
    by=np.fromfile(f4,dtype='float64').reshape(-1,2)
    bz=np.fromfile(f5,dtype='float64').reshape(-1,2)
    exsave[it,:]=ex[:,1];eysave[it,:]=ey[:,1];ezsave[it,:]=ez[:,1]
    bysave[it,:]=by[:,1];bzsave[it,:]=bz[:,1]

CPU times: user 641 ms, sys: 819 ms, total: 1.46 s
Wall time: 11.4 s
In [3]:
%%time
xstr=400
xend=1000
tmax=nt*400*0.05
fig=plt.figure(figsize=(15,6))
ax1 = fig.add_subplot(231)
ax2 = fig.add_subplot(232, sharex=ax1, sharey=ax1)
ax3 = fig.add_subplot(233, sharex=ax1, sharey=ax1)
ax4 = fig.add_subplot(235, sharex=ax1, sharey=ax1)
ax5 = fig.add_subplot(236, sharex=ax1, sharey=ax1)
im1=ax1.imshow(exsave,aspect='auto',origin='lower',cmap='bwr',extent=[0,1000,0,tmax],clim=(-0.2,0.2))
im2=ax2.imshow(eysave,aspect='auto',origin='lower',cmap='terrain_r',extent=[0,1000,0,tmax])
im3=ax3.imshow(ezsave,aspect='auto',origin='lower',cmap='bwr',extent=[0,1000,0,tmax],clim=(-0.01,0.01))
im4=ax4.imshow(bysave,aspect='auto',origin='lower',cmap='bwr',extent=[0,1000,0,tmax],clim=(-0.01,0.01))
im5=ax5.imshow(bzsave,aspect='auto',origin='lower',cmap='terrain_r',extent=[0,1000,0,tmax],clim=(0,3))
ax1.set_xlabel(r'$x/(c/\omega_{e})$');ax4.set_xlabel(r'$x/(c/\omega_{e})$');ax5.set_xlabel(r'$x/(c/\omega_{e})$')
ax1.set_ylabel(r'$t\omega_e$');ax4.set_ylabel(r'$t\omega_e$')
ax1.set_xlim(xstr,xend);
plt.colorbar(im1,ax=ax1)
plt.colorbar(im2,ax=ax2)
plt.colorbar(im3,ax=ax3)
plt.colorbar(im4,ax=ax4)
plt.colorbar(im5,ax=ax5)
plt.show()

CPU times: user 5.65 s, sys: 675 ms, total: 6.33 s
Wall time: 6.65 s
Shock reformation

We show the animation of a electron and ion phase space plot and the electric field component $E_x$.
At the foot region, there is a ion beam reflected from the shock front.
Due to the ion beam, electrons are heated.
Here is an animation of the shock reformation.

In [4]:
plt.rcParams['animation.html'] = 'html5'
In [5]:
ist=340
iet=410
ims=[]
fig=plt.figure(figsize=(10,5))
ax1=fig.add_subplot(111)
ax2=ax1.twinx()

for it in range(ist,iet,2):
    f1=open('xp1_%05d' %(it),'rb') 
    f2=open('xp2_%05d' %(it),'rb')
    f3=open('ex_%05d' %(it),'rb')
    xp1=np.fromfile(f1,dtype='float64').reshape(-1,4)
    xp2=np.fromfile(f2,dtype='float64').reshape(-1,4)
    bz =np.fromfile(f3,dtype='float64').reshape(-1,2)
    
    im1=ax1.plot(xp1[:,0],xp1[:,1],'c,')
    im2=ax1.plot(xp2[:,0],xp2[:,1],"r,")
    ax1.set_ylabel(r'$v_x/c$')
    ax1.set_ylim(-0.5,0.5)
    
    im3=ax2.plot(bz[:,0],bz[:,1],'k-',lw=1)
    ax2.set_ylabel(r'$E_x$')
    ax2.set_ylim(-0.8,0.8)

    plt.xlabel(r'$x/(c/\omega_{e})$')
    plt.xlim(550,700)

    ims.append(im1+im2+im3)
    #plt.show()
    
plt.close()
In [6]:
#fig=plt.figure()
ani = animation.ArtistAnimation(fig, ims)
ani
Out[6]:
Thermalization of electrons

The condition for the Bunemann instability is $$\frac{(4M_A/3)^2}{\beta_e}\frac{m_e}{m_i}>1.$$ With the present simulation parameters, $$\frac{(4M_A/3)^2}{\beta_e}\frac{m_e}{m_i}\sim 8>1.$$

Therefore, the mechanism of the electron heating is the Bunemann instaiblity.

Fluid values (density, velocity, and temperature)
In [7]:
it=350
bs=1000
f1=open('xp1_%05d' %(it),'rb') 
f2=open('xp2_%05d' %(it),'rb')
f3=open('bz_%05d' %(it),'rb')
xp1=np.fromfile(f1,dtype='float64').reshape(-1,4)
xp2=np.fromfile(f2,dtype='float64').reshape(-1,4)
bz =np.fromfile(f3,dtype='float64').reshape(-1,2)

(den1,bns)=np.histogram(xp1[:,0],bins=bs)
(den2,bns)=np.histogram(xp2[:,0],bins=bs)
(avl1,bns)=np.histogram(xp1[:,0],bins=bs,weights=xp1[:,1])
(avl2,bns)=np.histogram(xp2[:,0],bins=bs,weights=xp2[:,1])
tmp1=0.5*((xp1[:,1]-avl1[np.int_(xp1[:,0])]/den1[np.int_(xp1[:,0])])**2+xp1[:,2]**2+xp1[:,3]**2)
tmp2=0.5*((xp2[:,1]-avl2[np.int_(xp2[:,0])]/den2[np.int_(xp2[:,0])])**2+xp2[:,2]**2+xp2[:,3]**2)
(temp1,bns)=np.histogram(xp1[:,0],bins=bs,weights=tmp1)
(temp2,bns)=np.histogram(xp2[:,0],bins=bs,weights=tmp2)
binx=0.5*(bns[1:]+bns[:-1])
plt.plot(xp1[:,0],xp1[:,1],',')
plt.plot(xp2[:,0],xp2[:,1],',')
plt.plot(bz[:,0],bz[:,1])
plt.plot(binx,temp1/den1*10)
plt.plot(binx,temp2/den2*10,'k')
plt.xlim(500,800)
plt.ylim(-0.5,1.0)
plt.show()
In [8]:
%%time
nt=485
bs=1000
den1save=np.zeros((nt,bs))
den2save=np.zeros((nt,bs))
avl1save=np.zeros((nt,bs))
avl2save=np.zeros((nt,bs))
tmp1save=np.zeros((nt,bs))
tmp2save=np.zeros((nt,bs))
for it in range(nt):
    #print(it)
    f1=open('xp1_%05d' %(it),'rb') 
    f2=open('xp2_%05d' %(it),'rb')
    xp1=np.fromfile(f1,dtype='float64').reshape(-1,4)
    xp2=np.fromfile(f2,dtype='float64').reshape(-1,4)
    (den1,bns)=np.histogram(xp1[:,0],bins=bs)
    (den2,bns)=np.histogram(xp2[:,0],bins=bs)
    (avl1,bns)=np.histogram(xp1[:,0],bins=bs,weights=xp1[:,1])
    (avl2,bns)=np.histogram(xp2[:,0],bins=bs,weights=xp2[:,1])
    tmp1=0.5*((xp1[:,1]-avl1[np.int_(xp1[:,0])]/den1[np.int_(xp1[:,0])])**2+xp1[:,2]**2+xp1[:,3]**2)
    tmp2=0.5*((xp2[:,1]-avl2[np.int_(xp2[:,0])]/den2[np.int_(xp2[:,0])])**2+xp2[:,2]**2+xp2[:,3]**2)
    (temp1,bns)=np.histogram(xp1[:,0],bins=bs,weights=tmp1)    
    (temp2,bns)=np.histogram(xp2[:,0],bins=bs,weights=tmp2)
    binx=0.5*(bns[1:]+bns[:-1])
    den1save[it,:]=den1
    den2save[it,:]=den2
    avl1save[it,:]=avl1/den1
    avl2save[it,:]=avl2/den2
    tmp1save[it,:]=temp1/den1
    tmp2save[it,:]=temp2/den2 

CPU times: user 1min 32s, sys: 24.3 s, total: 1min 56s
Wall time: 4min 22s
In [9]:
%%time
plt.figure(figsize=(15,5))
plt.subplot(231);plt.imshow(den1save,origin='lower',aspect='auto',cmap='terrain_r')
#plt.ylim(340,410)
plt.colorbar()
plt.subplot(232);plt.imshow(avl1save,origin='lower',aspect='auto',cmap='terrain_r')
#plt.ylim(340,410)
plt.colorbar()
plt.subplot(233);plt.imshow(tmp1save,aspect='auto',origin='lower',cmap='terrain_r')
#plt.ylim(340,410)
plt.colorbar()
#plt.clim(-0.,2)
plt.subplot(234);plt.imshow(den2save,aspect='auto',origin='lower',cmap='terrain_r')
#plt.ylim(340,410)
plt.colorbar()
#plt.clim(-0.2,0.2)
plt.subplot(235);plt.imshow(avl2save,origin='lower',aspect='auto',cmap='terrain_r')
#plt.ylim(340,410)
plt.colorbar()
plt.subplot(236);plt.imshow(tmp2save*400,aspect='auto',origin='lower',cmap='terrain_r')
#plt.ylim(340,410)
plt.colorbar()
plt.show()

CPU times: user 2.54 s, sys: 24.3 ms, total: 2.57 s
Wall time: 2.91 s
Energy spectrum
In [10]:
it  =400
xmin=700
xmax=800
bins=200
enmin=1e-6
enmax=10
f1=open('xp1_%05d' %(it),'rb') 
f2=open('xp2_%05d' %(it),'rb')
xp1=np.fromfile(f1,dtype='float64').reshape(-1,4)
xp2=np.fromfile(f2,dtype='float64').reshape(-1,4)
plt.figure(figsize=(15,5))
plt.subplot(121)
plt.plot(xp1[:,0],xp1[:,1],',')
plt.plot(xp2[:,0],xp2[:,1],',')
plt.fill([xmin,xmax,xmax,xmin],[-1,-1,1,1],color="k",alpha=0.4)
plt.ylim(-1,1)

ene1=1./np.sqrt(1.-xp1[:,1]**2-xp1[:,2]**2-xp1[:,3]**2)-1.
ene2=1./np.sqrt(1.-xp2[:,1]**2-xp2[:,2]**2-xp2[:,3]**2)-1.
plt.subplot(122)
plt.hist(ene1[np.where((xp1[:,0]>=xmin)&(xp1[:,0]<xmax))],bins=np.logspace(np.log10(enmin),np.log10(enmax),bins),log=True,histtype='step',color='b')
#plt.hist(ene1[:],bins=np.logspace(np.log10(1e-6),np.log10(1.0),200),log=True,histtype='step',color='b')
plt.hist(ene2[np.where((xp2[:,0]>=xmin)&(xp2[:,0]<xmax))],bins=np.logspace(np.log10(enmin),np.log10(enmax),bins),log=True,histtype='step',color='orange')
#plt.hist(ene2[:],bins=np.logspace(np.log10(1e-6),np.log10(1.0),200),log=True,histtype='step',color='orange')
plt.yscale('log')
plt.xscale('log')
plt.xlim(enmin,enmax)
plt.show()
Trajectory of particles
In [11]:
%%time
nt=400
trj1=np.zeros((nt,4))
trj2=np.zeros((nt,4))
for it in range(nt):
    f1=open('xp1_%05d' %(it),'rb') 
    f2=open('xp2_%05d' %(it),'rb')
    trj1[it,:]=np.fromfile(f1,dtype='float64').reshape(-1,4)[40000]
    trj2[it,:]=np.fromfile(f2,dtype='float64').reshape(-1,4)[50000]

CPU times: user 1.16 s, sys: 12 s, total: 13.1 s
Wall time: 2min 11s
In [12]:
plt.imshow(bzsave,aspect='auto',origin='lower',cmap='terrain_r',extent=[0,1000,0,400])
tt=np.arange(400)
plt.plot(trj1[:,0],tt,'k')
plt.plot(trj2[:,0],tt,'k')
plt.colorbar()
plt.clim(0,2)
plt.xlim(400,1000)
plt.show()
Extra
In [13]:
it=80
bs=1000
f1=open('xp1_%05d' %(it),'rb') 
f2=open('xp2_%05d' %(it),'rb')
f3=open('bz_%05d' %(it),'rb')
xp1=np.fromfile(f1,dtype='float64').reshape(-1,4)
xp2=np.fromfile(f2,dtype='float64').reshape(-1,4)
bz =np.fromfile(f3,dtype='float64').reshape(-1,2)

(den1,bns)=np.histogram(xp1[:,0],bins=bs)
(den2,bns)=np.histogram(xp2[:,0],bins=bs)
(avl1,bns)=np.histogram(xp1[:,0],bins=bs,weights=xp1[:,1])
(avl2,bns)=np.histogram(xp2[:,0],bins=bs,weights=xp2[:,1])
tmp1=0.5*((xp1[:,1]-avl1[np.int_(xp1[:,0])]/den1[np.int_(xp1[:,0])])**2+xp1[:,2]**2+xp1[:,3]**2)
tmp2=0.5*((xp2[:,1]-avl2[np.int_(xp2[:,0])]/den2[np.int_(xp2[:,0])])**2+xp2[:,2]**2+xp2[:,3]**2)
(temp1,bns)=np.histogram(xp1[:,0],bins=bs,weights=tmp1)
(temp2,bns)=np.histogram(xp2[:,0],bins=bs,weights=tmp2)
binx=0.5*(bns[1:]+bns[:-1])
plt.plot(xp1[:,0],xp1[:,1],',')
plt.plot(xp2[:,0],xp2[:,1],',')
plt.plot(bz[:,0],bz[:,1])
plt.plot(binx,temp1/den1*10)
plt.plot(binx,temp2/den2*10,'k')
plt.plot(binx,den1/1000)
#plt.xlim(0,300)
plt.ylim(-0.5,1.0)
plt.show()