Gyro Surfing Acceleration

Reference: Kuramitsu2005PRL

$$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: $$\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].$$

Also, we add a electric potential modeled for a collisionless shock: $${E_x} = - a{\operatorname{sech} ^2}\left( {\frac{x}{l}} \right).$$

Simulation parameters in the shock rest frame: $k^\prime=0.1$, $c_p^\prime=4.5$, $v_u^\prime=5.5$, $v_t^\prime=1$, $b_w^\prime=0.5$, $2al=m_pv_u^2/2$, and $kl=1$.

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

In [16]:
%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 [17]:
it=400
f1 =open('data/up1%05d.d' %(it) ,'rb')
up1 =np.fromfile(f1, dtype='float64' ).reshape(-1,6)
plt.plot(up1[:,3],np.sqrt(up1[:,4]**2+up1[:,5]**2),',r')
plt.xlabel(r'$v_x/v_A$')
plt.ylabel(r'$v_\perp/v_A$')
plt.xlim(-8,8)
plt.ylim( 0,12)
plt.show()
In [33]:
it=0
f1 =open('data/up1%05d.d' %(it) ,'rb')
f2 =open('data/up2%05d.d' %(it) ,'rb')

up1 =np.fromfile(f1, dtype='float64' ).reshape(-1,6)
up2 =np.fromfile(f2, dtype='float64' ).reshape(-1,6)

iit=300
f  =open('data/up1%05d.d' %(iit) ,'rb')
upp=np.fromfile(f, dtype='float64' ).reshape(-1,6)

Psi0=np.arctan2(up1[:,4],up1[:,5])
mu0 =(up1[:,3]-5.5)/np.sqrt((up1[:,3]-5.5)**2+up1[:,4]**2+up1[:,5]**2)

plt.hist2d(Psi0,mu0,bins=(50,50),range=[[-np.pi,np.pi],[-1,1]],weights=np.sqrt(upp[:,4]**2+upp[:,5]**2),cmap='jet')
plt.colorbar()
plt.xlabel(r'$\Psi_0$')
plt.ylabel(r'$\mu_0$')
#plt.clim(0,50)
plt.show()