Cyclotron Resonance¶
Reference: Terasawa2012SpaceSci.Rev.
Equation of Motion¶
We solve the equation of motion for electrons using the Buneman-Boris method, $$m_e\frac{{d{\mathbf{v}}}}{{dt}} = -e\left( {{\mathbf{E}} + \frac{{\mathbf{v}}}{c} \times {\mathbf{B}}} \right).$$
Electromagnetic fields are defined as follows: \begin{gathered} {\mathbf{B}} = \left( {{B_0} + {B_{w,x}},{B_{w,y}},{B_{w,z}}} \right), \hfill \\ {\mathbf{E}} = \left( {{E_{w,x}},{E_{w,y}},{E_{w,z}}} \right). \hfill \\ \end{gathered} Here $B_0$ is the background magnetic field and EM fields with the subscript $w$ are wave components.
How to Determine the Wave Components¶
Deriving the dispersion realtion of a cold plasma, we get \begin{gathered} \tilde M{{\mathbf{E}}_w} = 0, \hfill \\ \tilde M = \left[ {\begin{array}{*{20}{c}} {P - {N^2}{{\sin }^2}\theta }&0&{{N^2}\sin \theta \cos \theta } \\ 0&{S - {N^2}}&{ - iD} \\ {{N^2}\sin \theta \cos \theta }&{iD}&{S - {N^2}{{\cos }^2}\theta } \end{array}} \right]. \hfill \\ \end{gathered}
Here, \begin{gathered} N = kc/\omega, \hfill \\ P = 1 - \frac{{\omega _{pi}^2}}{{{\omega ^2}}} - \frac{{\omega _{pe}^2}}{{{\omega ^2}}}, \hfill \\ R = 1 - \frac{{\omega _{pi}^2}}{{\omega \left( {\omega + {\Omega _{ci}}} \right)}} - \frac{{\omega _{pe}^2}}{{\omega \left( {\omega - {\Omega _{ce}}} \right)}}, \hfill \\ L = 1 - \frac{{\omega _{pi}^2}}{{\omega \left( {\omega - {\Omega _{ci}}} \right)}} - \frac{{\omega _{pe}^2}}{{\omega \left( {\omega + {\Omega _{ce}}} \right)}}, \hfill \\ S = \frac{1}{2}\left( {R + L} \right), \hfill \\ D = \frac{1}{2}\left( {R - L} \right). \hfill \\ \end{gathered}
The wavenumber and frequency $k$ and $\omega$ are determined to satisfy the dispersion relation.
Then we can calculate the amplitude of electric fields, \begin{gathered} {E_{w,x}} = - \frac{{{N^2}\sin \theta \cos \theta }}{{P - {N^2}{{\sin }^2}\theta }}{E_{w,z}}, \hfill \\ {E_{w,y}} = \frac{{iD}}{{S - {N^2}}}{E_{w,z}} \hfill \\ \end{gathered}
Using the definition below $${{\mathbf{E}}_w} \propto \exp i\left( {{\mathbf{k}} \cdot {\mathbf{r}} - \omega t} \right),$$ we can write down electric fields as follows $$\left( {{E_{w,x}},{E_{w,y}},{E_{w,z}}} \right) = \left( {{g_x}\cos \varphi ,{g_y}\sin \varphi ,{g_z}\cos \varphi } \right).$$
From the Faraday's law, we compute magnetic fields, \begin{gathered} {B_{w,x}} = - N{g_y}\sin \theta \sin \varphi , \hfill \\ {B_{w,y}} = N\left( {{g_x}\sin \theta - {g_z}\cos \theta } \right)\cos \varphi , \hfill \\ {B_{w,z}} = N{g_y}\cos \theta \sin \varphi . \hfill \\ \end{gathered}
Herafter, we pick up $Ng_y$ as the reference of the strength of waves.
def wave(th,k):
rm=1836
c=1e4
w=2.72e-3
b0=1e-4
N=k*c/w
wce=1
vA=1
wci=wce/rm
wpi=wce*c/vA/rm
wpe=wce*c/vA/np.sqrt(rm)
P=1-wpi**2/w**2-wpe**2/w**2
R=1-wpi**2/(w*(w+wci))-wpe**2/(w*(w-wce))
L=1-wpi**2/(w*(w-wci))-wpe**2/(w*(w+wce))
S=0.5*(R+L)
D=0.5*(R-L)
ex= N*np.sin(th)*np.cos(th)/(P-N**2*np.sin(th)**2)*(S-N**2)/D*b0
ey= b0/N#-D/(S-N**2)
ez=-(S-N**2)/D*b0/N
bx=-b0*np.sin(th)
by= N*(ex*np.cos(th)-ez*np.cos(th))
bz= b0*np.cos(th)
print(r'the amplitude of Ewx=',ex*c,'[B0(vA/c)]')
print(r'the amplitude of Ewy=',ey*c,'[B0(vA/c)]')
print(r'the amplitude of Ewz=',ez*c,'[B0(vA/c)]')
print(r'the amplitude of Bwx=',bx,'[B0]')
print(r'the amplitude of Bwy=',by,'[B0]')
print(r'the amplitude of Bwz=',bz,'[B0]')
%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
Simulation Results: Parallel Propagating Waves¶
Simulation Parameters:
- $\omega=5\Omega_{ci}=2.72\times 10^{-3}\Omega_{ce}$
- $k=2.04562\Omega_{ci}/v_A$
- $|B_w|/B_0=10^{-4}$
th=0 ;k=2.045618/5*2.72e-3
wave(th,k)
Particles are initially distributed as follows; $${\mathbf{r}} = \left( {0,0,0} \right),\,\,\,\,{\mathbf{v}} = \left( {{V_0},{V_ \bot }\cos \alpha ,{V_ \bot }\sin \alpha } \right),$$ where $V_0$ ranges from $-6000V_A$ to $6000V_A$ and $\alpha$ ranges from $0$ to $2\pi$.
We plot the variance of $\mu$, the cosine of the pitch angle, against the average velocity $<v_x>$. Definitions are described as follows; \begin{gathered} \left\langle {{v_x}} \right\rangle = \frac{1}{T}\int\limits_0^T {{v_x}\left( t \right)dt} , \hfill \\ \left\langle \mu \right\rangle = \frac{1}{T}\int\limits_0^T {\frac{{{v_x}\left( t \right)}}{{\left| {{\mathbf{v}}\left( t \right)} \right|}}dt} , \hfill \\ \left\langle {\delta {\mkern 1mu} \mu \delta \mu } \right\rangle = \left\langle \mu \right\rangle = \frac{1}{T}\int\limits_0^T {{{\left( {\frac{{{v_x}\left( t \right)}}{{\left| {{\mathbf{v}}\left( t \right)} \right|}} - \left\langle \mu \right\rangle } \right)}^2}dt} . \hfill \\ \end{gathered}
nt=4000; isav=1
nop=100000
dt=0.05*isav
%%time
vave=np.zeros(nop)
muave=np.zeros(nop)
for it in range(nt//isav):
f =open('data/up1%05d.d' %(it) ,'rb')
up =np.fromfile(f, dtype='float64' ).reshape(-1,6)
vave =vave+up[:,3]/(nt//isav)
mu =up[:,3]/np.sqrt(up[:,3]**2+up[:,4]**2+up[:,5]**2)
muave=muave+mu/(nt//isav)
mumuave=np.zeros(nop)
for it in range(nt//isav):
f =open('data/up1%05d.d' %(it) ,'rb')
up =np.fromfile(f, dtype='float64' ).reshape(-1,6)
mu =up[:,3]/np.sqrt(up[:,3]**2+up[:,4]**2+up[:,5]**2)
mumuave=mumuave+(mu-muave)**2/(nt//isav)
plt.plot(vave,np.sqrt(mumuave),',')
plt.yscale('log')
plt.xlabel(r'$<v_x>/V_A$')
plt.ylabel(r'$<\delta\mu\delta\mu>$')
plt.grid()
plt.show()
Simulation Results: Oblique Propagating Waves¶
Simulation Parameters:
- $\omega=5\Omega_{ci}=2.72\times 10^{-3}\Omega_{ce}$
- $k=2.42\Omega_{ci}/v_A$
- $|B_w|/B_0=10^{-4}$
nt=10000; isav=1
nop=50000
dt=0.05*isav
th=np.pi/4;k=2.42/5*2.72e-3
wave(th,k)
%%time
vave=np.zeros(nop)
muave=np.zeros(nop)
for it in range(nt//isav):
f =open('data2/up1%05d.d' %(it) ,'rb')
up =np.fromfile(f, dtype='float64' ).reshape(-1,6)
vave =vave+up[:,3]/(nt//isav)
mu =up[:,3]/np.sqrt(up[:,3]**2+up[:,4]**2+up[:,5]**2)
muave=muave+mu/(nt//isav)
mumuave=np.zeros(nop)
for it in range(nt//isav):
f =open('data2/up1%05d.d' %(it) ,'rb')
up =np.fromfile(f, dtype='float64' ).reshape(-1,6)
mu =up[:,3]/np.sqrt(up[:,3]**2+up[:,4]**2+up[:,5]**2)
mumuave=mumuave+(mu-muave)**2/(nt//isav)
plt.plot(vave,np.sqrt(mumuave),',')
plt.yscale('log')
plt.xlabel(r'$<v_x>/V_A$')
plt.ylabel(r'$<\delta\mu\delta\mu>$')
plt.grid()
plt.show()
