SPA: Heavy Ion Acceleration

Reference: Mizuta

$$m\frac{{d{\mathbf{v}}}}{{dt}} = q\left( {{\mathbf{E}} + \frac{{\mathbf{v}}}{c} \times {\mathbf{B}}} \right)$$

Normalization:

\begin{gathered} t \to t^\prime{\Omega _p} \hfill \\ x \to x^\prime \frac{{{v_A}}}{{{\Omega _p}}} \hfill \\ B,E \to B^\prime {B_0},E^\prime {B_0} \hfill \\ \end{gathered}

$$\frac{{d{\mathbf{v}}^\prime }}{{dt^\prime }} = \left( {\frac{c}{{{v_A}}}} \right){\mathbf{E}}^\prime + {\mathbf{v}}^\prime \times {\mathbf{B}}^\prime $$

Electromagnetic fields obey the Faraday's law:$${{\mathbf{k}}^\prime } \times {{\mathbf{E}}^\prime } = \left( {\frac{{{v_A}}}{c}} \right)\left( {\frac{\omega }{{{\Omega _p}}}} \right){{\mathbf{B}}^\prime }.$$

Electromagnetic fields are: $${\mathbf{B}} = \left( {{B_0},{B_y},{B_z}} \right),\,\,{\mathbf{E}} = \left( {0,{E_y},{E_z}} \right).$$

We put two circulary polarized waves (SPA and SBA): $$\left( {\begin{array}{*{20}{c}} {{B_y}} \\ {{B_z}} \end{array}} \right) = \sum\limits_{j = 1}^2 {{B_{ \bot ,j}}\left( {\begin{array}{*{20}{c}} {\sin \left( {k_j^\prime {x^\prime } - \omega _j^\prime {t^\prime } + {\phi _j}} \right)} \\ {\cos \left( {k_j^\prime {x^\prime } - \omega _j^\prime {t^\prime } + {\phi _j}} \right)} \end{array}} \right)}. $$

Simulation parameters in the shock rest frame: $B_{\perp1}^\prime=0.022$, $B_{\perp2}^\prime=0.1$, $V_{ph1}^\prime=18$, $V_{ph2}^\prime=-0.83$, $v_u^\prime=1$, $v_t^\prime=0.1$.

The wave number and frequency are determimed from the resonant condition: $$\left( {\omega + {\Omega _\alpha }} \right)/k = {V_R}.$$ Here $V_R=V_A$.

The parameter $v_u^\prime$ is the drift velocity.

We consider $\alpha$ particles.

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 IPython import display
plt.rcParams['font.size'] = 14
plt.rcParams['axes.linewidth'] = 1.5
In [2]:
nt=2000
isav=10
dt=0.05
In [3]:
it=200
f2 =open('data/up2%05d.d' %(it) ,'rb')
up2 =np.fromfile(f2, dtype='float64' ).reshape(-1,6)
In [7]:
th=np.linspace(0,2*np.pi,100)
for ir in range(25):
    plt.plot(ir*2*(np.cos(th))+18  ,ir*2*np.sin(th),'k',lw=0.5)
    plt.plot(ir*2*(np.cos(th))-0.83,ir*2*np.sin(th),'k',lw=0.5)

plt.plot(up2[:,3],np.sqrt(up2[:,4]**2+up2[:,5]**2),',')
plt.xlabel(r'$v_x/v_A$'); plt.ylabel(r'$v_\perp/v_A$')
plt.xlim(-25,25)
plt.ylim(  0,25)
plt.show()
In [14]:
nv=200
vp2hst =np.zeros((nt//isav,nv))
mu2hst =np.zeros((nt//isav,nv))
for it in range(nt//isav):
    f2 =open('data/up2%05d.d' %(it) ,'rb')
    up2 =np.fromfile(f2, dtype='float64' ).reshape(-1,6)
    uperp2=np.sqrt(up2[:,4]**2+up2[:,5]**2)
    mu2   =up2[:,3]/np.sqrt(up2[:,3]**2+up2[:,4]**2+up2[:,5]**2)
    vp2hst[it,:], bins=np.histogram(uperp2,bins=nv,range=[0,20])
    mu2hst[it,:], bins=np.histogram(mu2   ,bins=nv,range=[-1,1])
In [12]:
tt=np.arange(nt)*dt
plt.imshow(vp2hst.T,origin='lower',aspect='auto',cmap='jet',extent=[0,nt*dt,0,20])
plt.colorbar()
plt.clim(0,4000)
plt.xlabel(r'$t\Omega_p$'); plt.ylabel(r'$v_\perp/v_A$')
plt.show()
In [24]:
plt.imshow(mu2hst.T,origin='lower',aspect='auto',cmap='jet',extent=[0,nt*dt,0,20])
plt.colorbar()
plt.clim(0,5000)
plt.xlabel(r'$t\Omega_p$'); plt.ylabel(r'$\mu$')
plt.show()