Scatter Free Ion Acceleration

Reference: Sugiyama2001JGR

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

We put a circulary polarized wave in the background magnetic field $\bf{B}_0=B_0\bf{e}_x$: $$\tilde b = {b_y} + i{b_z} = {b_w}\exp \left[ {i\left( {k^\prime x^\prime - \omega^\prime t^\prime + {\phi _w}} \right)} \right].$$

This wave is compressed and amplified at the shock front and the shock profile is modeled as: $$1 + \frac{1}{2}\left( {r - 1} \right)\left[ {1 + \tanh \left( {\frac{x}{D}} \right)} \right].$$

Simulation parameters in the shock rest frame: $k^\prime=2\pi/56.6$, $c_p^\prime=4$, $v_u^\prime=5$, $v_t^\prime=1.6$, $b_w^\prime=0.6$.

We assume a shell distribution with the thermal velocity $v_t^\prime$.

%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

Model Description

lw=56.6
nx=1000
D=0.01*56.6
r=4
x=np.linspace(-50,50,nx)
ysh=1+0.5*(r-1)*(1+np.tanh(x/D))
kw=2*np.pi/lw*ysh
yw=np.sin(kw*x)*0.2*ysh+2
plt.plot(x,ysh,'-k')
plt.plot(x,yw)
plt.axis('off')
plt.show()

Simulation Results

it=1000
f1 =open('data/up1%05d.d' %(it) ,'rb')
up1 =np.fromfile(f1, dtype='float64' ).reshape(-1,6)

th=np.linspace(0,2*np.pi,100)
for ir in range(20):
    r=ir*1.0
    plt.plot(r*(np.cos(th))+4,r*np.sin(th),'k',lw=0.5)
    plt.plot(r*(np.cos(th))+1,r*np.sin(th),'--k',lw=0.5)
    
plt.plot(up1[:,3],np.sqrt(up1[:,4]**2+up1[:,5]**2),'.r',markersize=0.4)
plt.xlabel(r'$v_x/v_A$')
plt.ylabel(r'$v_\perp/v_A$')
plt.xlim(-12,12)
plt.ylim(  0,12)
plt.show()

Energy Spectrum

xmin=-100
xmax= 100
bins=200
enmin=1e1
enmax=1000

ene1=up1[:,3]**2+up1[:,4]**2+up1[:,5]**2
plt.hist(ene1[np.where((up1[:,0]>=xmin)&(up1[:,0]<xmax))],bins=np.logspace(np.log10(enmin),np.log10(enmax),bins),log=True,histtype='step',color='b')
plt.yscale('log')
plt.xscale('log')
plt.xlabel('Energy')
plt.ylabel('df(E)/dE')
plt.xlim(enmin,enmax)
plt.show()

Trajectory of an accelerated ion

it=1000
f1 =open('data/up1%05d.d' %(it), 'rb')
up1=np.fromfile(f1,dtype='float64').reshape(-1,6)

print(np.where(up1[:,3]**2+up1[:,4]**2+up1[:,5]>120))

(array([ 1801, 45126, 80622]),)

nt=1000; isav=1
ptrj=np.zeros((nt//isav,3))
ip=45126
for it in range(nt//isav):
    f1 =open('data/up1%05d.d' %(it) ,'rb')
    up1 =np.fromfile(f1, dtype='float64' ).reshape(-1,6)

    ptrj[it,0]=up1[ip,0]
    ptrj[it,1]=up1[ip,3]
    ptrj[it,2]=np.sqrt(up1[ip,4]**2+up1[ip,5]**2)

plt.plot(ptrj[:,0],ptrj[:,1]**2+ptrj[:,2]**2)
plt.xlabel(r'$x/(v_A/\Omega_p)$')
plt.ylabel('Energy')
plt.show()

th=np.linspace(0,2*np.pi,100)
for ir in range(20):
    r=ir*1.0
    plt.plot(r*(np.cos(th))+4,r*np.sin(th),'k',lw=0.5)
    plt.plot(r*(np.cos(th))+1,r*np.sin(th),'--k',lw=0.5)

plt.plot(ptrj[:,1],ptrj[:,2],'-r')
plt.xlim(-12,12)
plt.ylim(  0,12)
plt.xlabel(r'$v_x/v_A$')
plt.ylabel(r'$v_\perp/v_A$')
plt.show()