In [1]:
%config InlineBackend.figure_format = 'retina'
import numpy as np
import matplotlib.pyplot as plt
from IPython import display
import scipy.fftpack as sf
import matplotlib.animation as animation
import math as mt
import matplotlib as mpl
plt.rcParams['font.size'] = 14
plt.rcParams['axes.linewidth'] = 1.5
plt.rcParams['animation.html'] = 'jshtml'
plt.rcParams['animation.embed_limit'] = 2**128

Three Wave Coupling in a Hot Unmagnetized Plasma

D. A. Gurnett and A. Bhattacharjee, "Introduction to Plasma Physics (second edition)", Chap. 11

We consider a simplified 1D system described as follows, $$ \begin{align} \left(\partial_t^2-c^2\partial_x^2+\omega_{pe}^2\right)E_T &=-\omega_{pe}^2n_AE_L; \\ \left(\partial_t^2-C_L^2\partial_x^2+\omega_{pe}^2\right)E_L&=-\omega_{pe}^2n_AE_T; \\ \left(\partial_t^2-C_A^2\partial_x^2\right)n_A &=\frac{m_e}{m_i}\partial_x^2\left(U_TU_L\right); \\ \partial_tU_T&=-\frac{e}{m_e}E_T; \\ \partial_tU_L&=-\frac{e}{m_e}E_L. \end{align} $$

where, $E_T$ is the transverse electric field, $E_L$ the longitudinal electric field, $n_A$ the ion density fluctuation, and $U_{T(L)}$ the electron plasma velocity associated with $E_{T(L)}$. Here, $C_L^2=3\kappa T_e/m_e$ and $C_A^2=3\kappa T_e/m_i$. The right-hand side of the first three equations is a nonlinear term.

Linear Dispersion Relation: Transverse EM mode, Langmuir mode, and Ion Acoustic mode

Linearizing the first three equations, we obtain the following three linear dispersion relations, $$\omega^2=\omega_{pe}^2+c^2k^2 \quad \text{(Transverse EM mode)};$$ $$\omega^2=\omega_{pe}^2+C_L^2k^2 \quad \text{(Langmuir mode)};$$ $$\omega^2=C_A^2k^2 \quad \text{(Ion acoustic mode)}.$$

We plot the dispersion relations below with parameters, $m_i/m_e=100$, $C_L^2=0.05c^2$.
The time and length is normalized by $\omega_{pe}^{-1}$ and $c/\omega_{pe}$, respectively.

In [3]:
kk=np.linspace(-32,32,1000)
rm=100
CAe2=0.05
CAi2=CAe2/rm

wT=np.sqrt(1+kk**2)
wL=np.sqrt(1+CAe2*kk**2)
wA=np.sqrt(  CAi2*kk**2)

fig=plt.figure(figsize=(12,8))
plt.plot(kk,wT,label='EM')
plt.plot(kk,wL,label='Langmuir')
plt.plot(kk,wA,label='Ion Acoustic')
plt.legend()
plt.xlabel(r'$ck/\omega_{pe}$')
plt.ylabel(r'$\omega/\omega_{pe}$')
plt.xlim(-8,8)
plt.ylim(0,8)
plt.show()

Numerical Simulation

We solve the governing equations using the spectral method and the 4th order Runge-Kutta shceme.

In [4]:
def TWC(rm,nx,xmax,nt,dt,ET,ETt,EL,ELt,nA,nAt,uT,uL,isav,vte2,CA2):
    global kk,kkd

    ETf =sf.fft(ET )
    ETft=sf.fft(ETt)
    ELf =sf.fft(EL )
    ELft=sf.fft(ELt)
    nAf =sf.fft(nA )
    nAft=sf.fft(nAt)
    uTf =sf.fft(uT)
    uLf =sf.fft(uL)    

    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
 

    EThst=np.zeros((int(nt/isav),nx))
    ELhst=np.zeros((int(nt/isav),nx))
    nAhst=np.zeros((int(nt/isav),nx))
    uThst=np.zeros((int(nt/isav),nx))
    uLhst=np.zeros((int(nt/isav),nx))
    EThst[0,:]=np.real(sf.ifft(ETf))
    ELhst[0,:]=np.real(sf.ifft(ELf))
    nAhst[0,:]=np.real(sf.ifft(nAf))
    uThst[0,:]=np.real(sf.ifft(uTf))
    uLhst[0,:]=np.real(sf.ifft(uLf))
    for it in range(nt):
        
        ##---4th-order Runge-Kutta---#
        gt1,at1,gl1,al1,gA1,aA1,guT1,guL1=f(ETf           ,ETft           ,ELf           ,ELft           ,nAf           ,nAft           ,uTf           ,uLf           )
        gt2,at2,gl2,al2,gA2,aA2,guT2,guL2=f(ETf+0.5*dt*gt1,ETft+0.5*dt*at1,ELf+0.5*dt*gl1,ELft+0.5*dt*al1,nAf+0.5*dt*gA1,nAft+0.5*dt*aA1,uTf+0.5*dt*guT1,uLf+0.5*dt*guL1)
        gt3,at3,gl3,al3,gA3,aA3,guT3,guL3=f(ETf+0.5*dt*gt2,ETft+0.5*dt*at2,ELf+0.5*dt*gl2,ELft+0.5*dt*al2,nAf+0.5*dt*gA2,nAft+0.5*dt*aA2,uTf+0.5*dt*guT2,uLf+0.5*dt*guL2)
        gt4,at4,gl4,al4,gA4,aA4,guT4,guL4=f(ETf+    dt*gt3,ETft+    dt*at3,ELf+    dt*gl3,ELft+    dt*al3,nAf+    dt*gA3,nAft+    dt*aA3,uTf+    dt*guT3,uLf+    dt*guL3)
        ETf =ETf +dt*(gt1+2*gt2+2*gt3+gt4)/6
        ELf =ELf +dt*(gl1+2*gl2+2*gl3+gl4)/6
        nAf =nAf +dt*(gA1+2*gA2+2*gA3+gA4)/6
        uTf =uTf +dt*(guT1+2*guT2+2*guT3+guT4)/6
        uLf =uLf +dt*(guL1+2*guL2+2*guL3+guL4)/6
        ETft=ETft+dt*(at1+2*at2+2*at3+at4)/6
        ELft=ELft+dt*(al1+2*al2+2*al3+al4)/6
        nAft=nAft+dt*(aA1+2*aA2+2*aA3+aA4)/6
        
        if it%isav==0: 
            EThst[it//isav,:]=np.real(sf.ifft(ETf))
            ELhst[it//isav,:]=np.real(sf.ifft(ELf))
            nAhst[it//isav,:]=np.real(sf.ifft(nAf))
            uThst[it//isav,:]=np.real(sf.ifft(uTf))
            uLhst[it//isav,:]=np.real(sf.ifft(uLf))
            
    return locals()

def f(ETf,ETft,ELf,ELft,nAf,nAft,uTf,uLf):
    ### Linear terms ###
    aETf =ETft
    aETft=(1j*kk)**2*ETf-ETf
    aELf =ELft
    aELft=vte2*(1j*kk)**2*ELf-ELf
    anAf =nAft
    anAft=CA2*(1j*kk)**2*nAf
    auTf =-ETf
    auLf =-ELf
    
    #### nonlinear terms ###
    ET=np.real(sf.ifft(ETf*kkd))
    EL=np.real(sf.ifft(ELf*kkd))
    nA=np.real(sf.ifft(nAf*kkd))
    uT=np.real(sf.ifft(uTf*kkd))
    uL=np.real(sf.ifft(uLf*kkd))
    aETft=aETft+sf.fft(-nA*EL)
    aELft=aELft+sf.fft(-nA*ET)
    anAft=anAft+(1j*kk)**2*sf.fft(uT*uL)/rm
 
    return aETf,aETft,aELf,aELft,anAf,anAft,auTf,auLf

Linear modes

First we check linear modes in this system by setting small amplitude waves (or fluctuations).
The parameters are $N_x=64$, $\delta x=2\pi/N_x$, $\Delta t=0.1*\Delta x$, $m_i/m_e=100$, and $C_L^2=0.05$.

In [5]:
nx=64
xmax=2*np.pi
dx=xmax/nx
dt=0.1*dx
isav=2
nt=8192*2
rm=100
vte2=0.05
CA2 =vte2/rm

#---initial condition---#
xx=dx*np.arange(nx)
ET =np.zeros(nx)
ETt=np.zeros(nx)
EL =np.zeros(nx)
ELt=np.zeros(nx)
nA =np.zeros(nx)
nAt=np.zeros(nx)
uT =np.zeros(nx)
uL =np.zeros(nx)

for ik in range(-32,32):
    bw  =1e-4*np.random.randn()
    phiw=2*np.pi*np.random.rand(3)
    ET =ET +bw*np.cos(ik*2*np.pi*xx/xmax+phiw[0])#+1j*bw*np.sin(ik*2*np.pi*xx/xmax+phiw)
    #ETt=ETt#-bw*np.sin(ik*2*np.pi*xx/xmax+phiw)#+1j*bw*np.sin(ik*2*np.pi*xx/xmax+phiw)
    EL =EL +bw*np.cos(ik*2*np.pi*xx/xmax+phiw[1])#+1j*bw*np.sin(ik*2*np.pi*xx/xmax+phiw)
    #ELt=ELt#-bw*np.sin(ik*2*np.pi*xx/xmax+phiw)#+1j*bw*np.sin(ik*2*np.pi*xx/xmax+phiw)
    nA =nA +bw*np.cos(ik*2*np.pi*xx/xmax+phiw[2])#+1j*bw*np.sin(ik*2*np.pi*xx/xmax+phiw)
    #nAt=nAt#-bw*np.sin(ik*2*np.pi*xx/xmax+phiw)#+1j*bw*np.sin(ik*2*np.pi*xx/xmax+phiw)
    
data=TWC(rm,nx,xmax,nt,dt,ET,ETt,EL,ELt,nA,nAt,uT,uL,isav,vte2,CA2)
locals().update(data)
In [6]:
kmin=2*mt.pi/(dx*nx)*(-nx/2)
kmax=2*mt.pi/(dx*nx)*(nx/2)
wmin=2*mt.pi/(dt*nt*isav)*(-nt/2)
wmax=2*mt.pi/(dt*nt*isav)*(nt/2)
kaxis=np.linspace(kmin,kmax,nx)
waxis=np.linspace(wmin,wmax,nt)

fftET=sf.fftshift(sf.fft2(EThst[:,::-1]))
fftEL=sf.fftshift(sf.fft2(ELhst[:,::-1]))
fftnA=sf.fftshift(sf.fft2(nAhst[:,::-1]))
plt.imshow(abs(fftET)+abs(fftEL)+abs(fftnA),origin='lower',aspect='auto',extent=([kmin,kmax,wmin,wmax]))
plt.axis([kmin,kmax,0,wmax/16])
plt.clim(0,8e0)
plt.colorbar()

kk=np.linspace(-32,32,1000)
wT=np.sqrt(1+     kk**2)
wL=np.sqrt(1+vte2*kk**2)
wA=np.sqrt(   CA2*kk**2)

plt.plot(kk,wT,'w--',lw=1)
plt.plot(kk,wL,'w--',lw=1)
plt.plot(kk,wA,'w--',lw=1)
plt.xlabel(r'$ck/\omega_{pe}$')
plt.ylabel(r'$\omega/\omega_{pe}$')
plt.show()
#plt.subplot(143)

Parametric Instability

Next, we inject a large-amplitude EM wave and observe how it decays into Langmuir and ion acoustic modes.
The parameters are $N_x=256$, $\Delta x=2\pi/N_x$, $\Delta t=0.1*\Delta x$, $m_i/m_e=100$, and $C_L^2=0.05$. The wave amplitude of the EM wave is $0.8$ and the wave number is $2$.

In [7]:
nx=256
xmax=2*np.pi
dx=xmax/nx
dt=0.04*dx
isav=16
nt=8192*16
rm=100
vte2=0.05
CA2 =vte2/rm

#---initial condition---#
xx=dx*np.arange(nx)
b0=0.8
phiw=2*np.pi*np.random.rand()
ET =np.zeros(nx)+                b0*np.cos(2*2*np.pi*xx/xmax+phiw)
ETt=np.zeros(nx)+np.sqrt(1+2**2)*b0*np.sin(2*2*np.pi*xx/xmax+phiw)
EL =np.zeros(nx)
ELt=np.zeros(nx)
nA =np.zeros(nx)
nAt=np.zeros(nx)
uT =np.zeros(nx)
uL =np.zeros(nx)

for ik in range(-128,128):
    bw  =1e-4*np.random.randn()
    phiw=2*np.pi*np.random.rand(3)
    ET =ET +bw*np.cos(ik*2*np.pi*xx/xmax+phiw[0])#+1j*bw*np.sin(ik*2*np.pi*xx/xmax+phiw)
    #ETt=ETt#-bw*np.sin(ik*2*np.pi*xx/xmax+phiw)#+1j*bw*np.sin(ik*2*np.pi*xx/xmax+phiw)
    EL =EL +bw*np.cos(ik*2*np.pi*xx/xmax+phiw[1])#+1j*bw*np.sin(ik*2*np.pi*xx/xmax+phiw)
    #ELt=ELt#-bw*np.sin(ik*2*np.pi*xx/xmax+phiw)#+1j*bw*np.sin(ik*2*np.pi*xx/xmax+phiw)
    nA =nA +bw*np.cos(ik*2*np.pi*xx/xmax+phiw[2])#+1j*bw*np.sin(ik*2*np.pi*xx/xmax+phiw)
    #nAt=nAt#-bw*np.sin(ik*2*np.pi*xx/xmax+phiw)#+1j*bw*np.sin(ik*2*np.pi*xx/xmax+phiw)
    
data=TWC(rm,nx,xmax,nt,dt,ET,ETt,EL,ELt,nA,nAt,uT,uL,isav,vte2,CA2)
locals().update(data)
In [8]:
def update_anim(it):
    
    fig.clf()

    ax1 = fig.add_subplot(111)
    
    ax1.clear() 
    
    xx=dx*np.arange(nx)
    im1=ax1.plot(xx,EThst[it*32,:],label=r'$E_T$')
    im2=ax1.plot(xx,ELhst[it*32,:],label=r'$E_L$')
    im3=ax1.plot(xx,nAhst[it*32,:],label=r'$n_A$')
    ax1.legend(loc='upper right')
    #ax1.set_title(r'$t=$%d' %(it*tsav))
    ax1.set_xlabel(r'$X/(c/\omega_{pe})$')
    ax1.set_ylabel(r'$\Psi$')
    ax1.set_xlim(0,nx*dx)
    ax1.set_ylim(-2,2)

fig=plt.figure(figsize=(16,4))
anim=animation.FuncAnimation(fig,update_anim,frames=nt//isav//32)
plt.close()
anim
Out[8]:
In [9]:
fig, ax = plt.subplots(ncols=3, sharex=True, figsize=(16, 4))
ax[0].imshow(EThst,cmap='bwr',origin='lower',aspect='auto',extent=([0,nx*dx,0,nt*dt]))
ax[1].imshow(ELhst,cmap='bwr',origin='lower',aspect='auto',extent=([0,nx*dx,0,nt*dt]))
ax[2].imshow(nAhst,cmap='bwr',origin='lower',aspect='auto',extent=([0,nx*dx,0,nt*dt]))
ax[0].set_title(r'$E_T$')
ax[1].set_title(r'$E_L$')
ax[2].set_title(r'$n_A$')
for i in range(3):
    ax[i].set_xlabel(r'$X/(c/\omega_{pe})$')
ax[0].set_ylabel(r'$T\omega_{pe}$')
plt.show()

The decay process occurs where the matching condition ($k_T=k_L+k_A$ and $\omega_T=\omega_L+\omega_A$) is satisfied.

In [10]:
kmin=2*mt.pi/(dx*nx)*(-nx/2)
kmax=2*mt.pi/(dx*nx)*(nx/2)
wmin=2*mt.pi/(dt*nt*isav)*(-nt/2)
wmax=2*mt.pi/(dt*nt*isav)*(nt/2)
kaxis=np.linspace(kmin,kmax,nx)
waxis=np.linspace(wmin,wmax,nt)

fftET=sf.fftshift(sf.fft2(EThst[:,::-1]))
fftEL=sf.fftshift(sf.fft2(ELhst[:,::-1]))
fftnA=sf.fftshift(sf.fft2(nAhst[:,::-1]))
plt.figure(figsize=(8,6))
plt.imshow(abs(fftET)+abs(fftEL)+abs(fftnA),cmap='gray_r',origin='lower',aspect='auto',extent=([kmin,kmax,wmin,wmax]))
plt.axis([kmin/8,kmax/8,0,wmax/64])
plt.clim(0,4e5)
plt.colorbar()
kk=np.linspace(-32,32,1000)
wT=np.sqrt(1+kk**2)
wL=np.sqrt(1+vte2*kk**2)
wA=np.sqrt(  CA2*kk**2)

plt.plot(kk,wT,'k--',lw=1)
plt.plot(kk,wL,'k--',lw=1)
plt.plot(kk,wA,'k--',lw=1)
plt.xlabel(r'$ck/\omega_{pe}$')
plt.ylabel(r'$\omega/\omega_{pe}$')

plt.arrow(0, 0, 2, np.sqrt(1+2**2), head_width=0.01, head_length=0.01, fc='C1', ec='C1')
plt.arrow(0, 0, 8.2,np.sqrt(1+vte2*8.2**2),head_width=0.01, head_length=0.01, fc='C1', ec='C1')
plt.arrow(0, 0,-6,np.sqrt(CA2*6**2),head_width=0.01, head_length=0.01, fc='C1', ec='C1')
plt.arrow(2,np.sqrt(1+2**2), 6,-np.sqrt(CA2*6**2),head_width=0.01, head_length=0.01, fc='C1', ec='C1',ls='--',width=0.001)
plt.arrow(2,np.sqrt(1+2**2),-8.2,-np.sqrt(1+vte2*8.2**2),head_width=0.01, head_length=0.01, fc='C1', ec='C1',ls='--',width=0.001)

plt.arrow(0, 0, 2, np.sqrt(1+2**2), head_width=0.1, head_length=0.5,length_includes_head=True, fc='C2', ec='C2',width=0.01)
plt.arrow(0, 0,-7.8,np.sqrt(1+vte2*7.8**2),head_width=0.1, head_length=0.5,length_includes_head=True, fc='C2', ec='C2',width=0.01)
plt.arrow(0, 0,9.8,np.sqrt(CA2*9.8**2),head_width=0.01, head_length=0.1,length_includes_head=True, fc='C2', ec='C2')
plt.arrow(2,np.sqrt(1+2**2),-9.8,-np.sqrt(CA2*9.8**2),head_width=0.01, head_length=0.01, fc='C2', ec='C2',ls='--',width=0.001)
plt.arrow(2,np.sqrt(1+2**2),7.8,-np.sqrt(1+vte2*7.8**2),head_width=0.01, head_length=0.01, fc='C2', ec='C2',ls='--',width=0.001)

plt.show()
In [2]:
!jupyter nbconvert --to html --template tmp.tpl Three_Wave_Coupling\(EM_Langmuir_Acoustic\).ipynb

[NbConvertApp] Converting notebook Three_Wave_Coupling(EM_Langmuir_Acoustic).ipynb to html
[NbConvertApp] Writing 24262590 bytes to Three_Wave_Coupling(EM_Langmuir_Acoustic).html