Nonlinear Evolution of Multi-ion Magnetosonic Waves

Reference: [Toida1995PoP]

We simulate the nonlinear evolution of magnetosonic waves in a multi-ion species plasma. The governing equations are described as follows, and we consider electron, proton, and alpha particle,

\begin{align} &{\partial _t}{n_j} + \nabla \cdot \left( {{n_j}{{\bf{v}}_j}} \right) = 0;\\ &{{\partial _t}{{\bf{v}}_j} + \left( {{{\bf{v}}_j} \cdot \nabla } \right)}{\bf{v}}_j = \frac{q_j}{m_j}\left( {{\bf{E}} + \frac{{{{\bf{v}}_j}}}{c} \times {\bf{B}}} \right);\\ &\frac{1}{c}{\partial _t}{\bf{B}} = - \nabla \times {\bf{E}};\\ &\frac{1}{c}{\partial _t}{\bf{E}} = \nabla \times {\bf{B}} - \frac{{4\pi }}{c}\sum\limits_j {{q_j}{n_j}{{\bf{v}}_j}}. \end{align}

Linear Dispersion Relation

Neglecting the displacement current, we obtain the dispersion relation for perpendicular propagating waves, $$ \left(\sum_j\frac{\Omega_{cj}\omega_{pj}^2}{\omega^2-\Omega_{cj}^2}\right)^2-\left(c^2k^2+\sum_j\frac{\omega_{pj}^2\omega^2}{\omega^2-\Omega_{cj}^2}\right)\left(\sum_j\frac{\omega_{pj}^2}{\omega^2-\Omega_{cj}^2}\right)=0 $$

In [1]:
%config InlineBackend.figure_format='retina'
import numpy as np
import matplotlib.pyplot as plt
import matplotlib as mpl
import scipy.fftpack as sf
import scipy as sp
import math as mt
import matplotlib.animation as animation
plt.rcParams['font.size'] = 12
plt.rcParams['axes.linewidth'] = 1.5
plt.rcParams['animation.embed_limit']=60

plt.rcParams['animation.html'] = 'jshtml'
In [2]:
q=np.zeros(3); m=np.zeros(3)

q[0]=-1; q[1]=1   ; q[2]=2
m[0]= 1; m[1]=1000; m[2]=4*m[1]

c=1
b0=0.5
ne0=1
na0=0.84
nb0=0.08

wpe=np.sqrt(ne0*q[0]**2/m[0])
wpa=np.sqrt(na0*q[1]**2/m[1])
wpb=np.sqrt(nb0*q[2]**2/m[2])
wce=q[0]*b0/m[0]
wca=q[1]*b0/m[1]
wcb=q[2]*b0/m[2]

wpe2=wpe**2
wpa2=wpa**2
wpb2=wpb**2
wce2=wce**2
wca2=wca**2
wcb2=wcb**2

kk=np.logspace(-3,1,1000)

A2=wpe2**2*(1+c**2*kk**2/wpe2)
A1=(wpa2/wca2+wpb2/wcb2)**2*wca2*wcb2*wce2+(wpa2+wpb2)*wce2*c**2*kk**2
A0=(wpa2/wca2+wpb2/wcb2)*wca2*wcb2*wce2*c**2*kk**2

wp2=(A1+np.sqrt(A1**2-4*A2*A0))/(2*A2)
wm2=(A1-np.sqrt(A1**2-4*A2*A0))/(2*A2)

wp=np.sqrt(wp2)
wm=np.sqrt(wm2)
plt.plot(kk,wp/wca,label=r'$\omega_{+}$')
plt.plot(kk,wm/wca,label=r'$\omega_{-}$')
plt.legend()
plt.xscale('log')
plt.yscale('log')
plt.xlim(1e-3,1e1)
plt.xlabel(r'$ck/\omega_{pe}$')
plt.ylabel(r'$\omega/\Omega_a$')
plt.show()

Nonlinear Simulation

We solve the governing equations using the pseudo-spectral method and 4th-order Runge-Kutta scheme.
We initially put a soliton solution with a large-amplitude $\delta B_0=0.2$ and impose a viscosity term for $v_{ex}$ to prevent oscillations in front of the shock.

In [8]:
def MIP(nx,nt,dx,dt,q,m,isav,n,vx,vy,ex,ey,bz):
    global kk,kkd
    
    xmax=nx*dx
    kk =2*np.pi/xmax*np.r_[np.arange(nx/2),np.arange(-nx/2,0)]
    kkd=np.zeros(nx)
    kkd[np.abs(kk)<=2/3*np.max(abs(kk))]=1.0 ### 2/3-rule for dealiasing
    
    nhst  =np.zeros((3,nt//isav+1,nx))
    vxhst =np.zeros((3,nt//isav+1,nx))
    vyhst =np.zeros((3,nt//isav+1,nx))
    exhst =np.zeros((nt//isav+1,nx))
    eyhst =np.zeros((nt//isav+1,nx))
    bzhst =np.zeros((nt//isav+1,nx))
    
    nhst[0,0,:]=n[0,:]
    nhst[1,0,:]=n[1,:]
    nhst[2,0,:]=n[2,:]
    vxhst[0,0,:]=vx[0,:]
    vxhst[1,0,:]=vx[1,:]
    vxhst[2,0,:]=vx[2,:]
    vyhst[0,0,:]=vy[0,:]
    vyhst[1,0,:]=vy[1,:]
    vyhst[2,0,:]=vy[2,:]
    exhst[0,:]=ex
    eyhst[0,:]=ey
    bzhst[0,:]=bz
    
    nf =sf.fft(n ,axis=1)
    vxf=sf.fft(vx,axis=1)
    vyf=sf.fft(vy,axis=1)
    exf=sf.fft(ex)
    eyf=sf.fft(ey)
    bzf=sf.fft(bz)
    
    for it in range(1,nt+1):
        ga1,gb1,gc1,gd1,ge1,gf1=f(nf           ,vxf           ,vyf           ,exf           ,eyf           ,bzf           )
        ga2,gb2,gc2,gd2,ge2,gf2=f(nf+0.5*dt*ga1,vxf+0.5*dt*gb1,vyf+0.5*dt*gc1,exf+0.5*dt*gd1,eyf+0.5*dt*ge1,bzf+0.5*dt*gf1)
        ga3,gb3,gc3,gd3,ge3,gf3=f(nf+0.5*dt*ga2,vxf+0.5*dt*gb2,vyf+0.5*dt*gc2,exf+0.5*dt*gd2,eyf+0.5*dt*ge2,bzf+0.5*dt*gf2)
        ga4,gb4,gc4,gd4,ge4,gf4=f(nf+    dt*ga3,vxf+    dt*gb3,vyf+    dt*gc3,exf+    dt*gd3,eyf+    dt*ge3,bzf+    dt*gf3)
        
        nf  =nf +dt*(ga1+2*ga2+2*ga3+ga4)/6
        vxf =vxf+dt*(gb1+2*gb2+2*gb3+gb4)/6
        vyf =vyf+dt*(gc1+2*gc2+2*gc3+gc4)/6
        exf =exf+dt*(gd1+2*gd2+2*gd3+gd4)/6
        eyf =eyf+dt*(ge1+2*ge2+2*ge3+ge4)/6
        bzf =bzf+dt*(gf1+2*gf2+2*gf3+gf4)/6
        
        #---add viscocity and resistivity term---#
        mu =1e-2*(1+(4*kk)**2)
        #for i in range(3):
        vxf[0,:]=np.exp(mu*(1j*kk)**2*dt)*vxf[0,:]
        #exf=np.exp(nu*(1j*kk)**2*dt)*exf
        
        if(it%isav==0): 
            nhst[ 0,it//isav,:]=np.real(sf.ifft(nf[ 0,:]))
            nhst[ 1,it//isav,:]=np.real(sf.ifft(nf[ 1,:]))
            nhst[ 2,it//isav,:]=np.real(sf.ifft(nf[ 2,:]))
            vxhst[0,it//isav,:]=np.real(sf.ifft(vxf[0,:]))
            vxhst[1,it//isav,:]=np.real(sf.ifft(vxf[1,:]))
            vxhst[2,it//isav,:]=np.real(sf.ifft(vxf[2,:]))
            vyhst[0,it//isav,:]=np.real(sf.ifft(vyf[0,:]))
            vyhst[1,it//isav,:]=np.real(sf.ifft(vyf[1,:]))
            vyhst[2,it//isav,:]=np.real(sf.ifft(vyf[2,:]))
            exhst[  it//isav,:]=np.real(sf.ifft(exf    ))
            eyhst[  it//isav,:]=np.real(sf.ifft(eyf    ))
            bzhst[  it//isav,:]=np.real(sf.ifft(bzf    ))

    return locals()


def f(nf,vxf,vyf,exf,eyf,bzf):
    an =np.zeros((3,nx),dtype='complex')
    avx=np.zeros((3,nx),dtype='complex')
    avy=np.zeros((3,nx),dtype='complex')
    aex=np.zeros(nx,dtype='complex')
    aey=np.zeros(nx,dtype='complex')
    abz=np.zeros(nx,dtype='complex')
    
    ### linear terms ###
    for i in range(3):
        avx[i,:]=q[i]/m[i]*exf[:]
        avy[i,:]=q[i]/m[i]*eyf[:]
    abz=-1j*kk*eyf
    aey=-1j*kk*bzf
    
    ### non-linear terms ###    
    n =np.real(sf.ifft(nf *kkd,axis=1))
    vx=np.real(sf.ifft(vxf*kkd,axis=1))
    vy=np.real(sf.ifft(vyf*kkd,axis=1))
    ex=np.real(sf.ifft(exf*kkd))
    ey=np.real(sf.ifft(eyf*kkd))
    bz=np.real(sf.ifft(bzf*kkd))
    vxx=np.real(sf.ifft(1j*kk*vxf*kkd,axis=1))
    vyx=np.real(sf.ifft(1j*kk*vyf*kkd,axis=1))
    
    jx=np.zeros(nx)
    jy=np.zeros(nx)
    for i in range(3):
        an[i,:] =-1j*kk*sf.fft(n[i,:]*vx[i,:])
        avx[i,:]=avx[i,:]-sf.fft(vx[i,:]*vxx[i,:])
        avy[i,:]=avy[i,:]-sf.fft(vx[i,:]*vyx[i,:])
        avx[i,:]=avx[i,:]+q[i]/m[i]*sf.fft(vy[i,:]*bz)
        avy[i,:]=avy[i,:]-q[i]/m[i]*sf.fft(vx[i,:]*bz)
        jx=jx+q[i]*n[i,:]*vx[i,:]
        jy=jy+q[i]*n[i,:]*vy[i,:]
    
    aex=sf.fft(-jx)
    aey=aey+sf.fft(-jy)
    
    return an,avx,avy,aex,aey,abz
In [9]:
%%time
nx=2048;nt=60000;dx=2;dt=1;isav=600
q=np.zeros(3); m=np.zeros(3)

q[0]=-1; q[1]=1   ; q[2]=2
m[0]= 1; m[1]=1000; m[2]=4*m[1]

b0=0.5
ne0=1
na0=0.84
nb0=0.08

db=0.2

wce=q[0]*b0/m[0]
wca=q[1]*b0/m[1]
wcb=q[2]*b0/m[2]

wpe=1
wpa=np.sqrt(na0*q[1]**2/m[1])
wpb=np.sqrt(nb0*q[2]**2/m[2])

wpe2=wpe**2
wpa2=wpa**2
wpb2=wpb**2
wca2=wca**2
wcb2=wcb**2

wmr=np.sqrt((wpa2*wcb2+wpb2*wca2)/(wpa2+wpb2))
wp0=(wpa2/wca2+wpb2/wcb2)*np.sqrt(wca2*wcb2)*abs(wce)

ro0=na0*m[1]+nb0*m[2]
vA =b0/np.sqrt(ro0)


kc  =wmr/vA
rab =1-(wmr/wp0)**2
dl  =np.sqrt(rab)/kc
vp0 =vA
lmbd=dl
alp =1

x=np.arange(nx)*dx-nx*dx/4
vp=vp0*(1+alp*db/2)
D=2*lmbd/np.sqrt(alp*db)

b1  =db/np.cosh(x/D)**2
ey1 =b1#*vA*b0
vex1=b1#*vA
vax1=b1#*vA
vbx1=b1#*vA
ne1 =b1#*ne0
na1 =b1#*na0
nb1 =b1#*nb0

jj=np.arange(nx)
jp=np.r_[jj[0:nx-1]+1,0]
jm=np.r_[nx-1,jj[1:nx]-1]
#dxb1=(b1[jp]-b1[jm])/(2*dx)*dl

xmax=dx*nx
kk =2*np.pi/xmax*np.r_[np.arange(nx/2),np.arange(-nx/2,0)]
dxb1=np.real(sf.ifft(1j*kk*sf.fft(b1)))*dl

ex1=-(vA)**2*wmr/np.sqrt(rab)*(wpe2/wce**3+wpa2/wca**3+wpb2/wcb**3)*dxb1#*vA*b0
vey1=-(wmr/wce/np.sqrt(rab)*dxb1+ex1)#*vA
vay1=-(wmr/wca/np.sqrt(rab)*dxb1+ex1)#*vA
vby1=-(wmr/wcb/np.sqrt(rab)*dxb1+ex1)#*vA


n =np.zeros((3,nx))
vx=np.zeros((3,nx))
vy=np.zeros((3,nx))
ex=np.zeros(nx)
ey=np.zeros(nx)
bz=np.zeros(nx)

n[0,:]=ne0*(1+ne1)
n[1,:]=na0*(1+na1)
n[2,:]=nb0*(1+nb1)

vx[0,:]=vex1*vA
vx[1,:]=vax1*vA
vx[2,:]=vbx1*vA
vy[0,:]=vey1*vA
vy[1,:]=vay1*vA
vy[2,:]=vby1*vA

ex=ex1*vA*b0
ey=ey1*vA*b0
bz=b0*(1+b1)

data=MIP(nx,nt,dx,dt,q,m,isav,n,vx,vy,ex,ey,bz)
locals().update(data)

CPU times: user 3min 36s, sys: 566 ms, total: 3min 37s
Wall time: 3min 40s
In [10]:
def update_anim(it):
    
    fig.clf()

    ax=fig.subplots(6,1,sharex=True) #add subplots
    #ax1 = fig.add_subplot(111)
    
    #fig.tight_layout() #reduce spacing around the figure
    
    #ax1.clear()

    xx=dx*np.arange(nx)
    ax[0].plot(xx,nhst[0, it,:],label=r'$n_e$')
    ax[0].plot(xx,nhst[1, it,:],label=r'$n_a$')
    ax[0].plot(xx,nhst[2, it,:]*8,label=r'$8\times n_b$')
    ax[1].plot(xx,vxhst[0,it,:],label=r'$v_{ex}$')
    ax[1].plot(xx,vxhst[1,it,:],label=r'$v_{ax}$')
    ax[1].plot(xx,vxhst[2,it,:],label=r'$v_{bx}$')
    ax[2].plot(xx,vyhst[0,it,:],label=r'$v_{ey}$')
    ax[2].plot(xx,vyhst[1,it,:],label=r'$v_{ay}$')
    ax[2].plot(xx,vyhst[2,it,:],label=r'$v_{by}$')    
    ax[3].plot(xx,exhst[  it,:])
    ax[4].plot(xx,eyhst[  it,:])
    ax[5].plot(xx,bzhst[  it,:])
    ax[0].legend(loc='upper left')
    ax[1].legend(loc='upper left')
    ax[2].legend(loc='upper left')
    ax[5].set_xlabel(r'$X/(c/\omega_{pe})$')
    ax[0].set_ylabel(r'$n$')
    ax[1].set_ylabel(r'$v_x$')
    ax[2].set_ylabel(r'$v_y$')
    ax[3].set_ylabel(r'$E_x$')
    ax[4].set_ylabel(r'$E_y$')
    ax[5].set_ylabel(r'$B_z$')
    ax[5].set_xlim(500,2500)
    ax[0].set_ylim(np.min(nhst*8),np.max(nhst))
    ax[1].set_ylim(np.min(vxhst),np.max(vxhst))
    ax[2].set_ylim(np.min(vyhst),np.max(vyhst))
    ax[3].set_ylim(np.min(exhst),np.max(exhst))
    ax[4].set_ylim(np.min(eyhst),np.max(eyhst))
    ax[5].set_ylim(np.min(bzhst),np.max(bzhst))
    ax[0].set_title(r'$t=$%d$\omega_{pe}^{-1}$' %(it*isav*dt))
In [11]:
fig=plt.figure(figsize=(18,12))
anim=animation.FuncAnimation(fig,update_anim,frames=nt//isav+1)
plt.close()
anim
Out[11]:
In [7]:
!jupyter nbconvert --to html --template tmp.tpl multi-ion_fluid.ipynb

[NbConvertApp] Converting notebook multi-ion_fluid.ipynb to html
[NbConvertApp] Writing 9884914 bytes to multi-ion_fluid.html