Kappa Distribution in Particle in Cell Simulation

reference: Zenitani & Nakano

Introduction

Kappa distribution is described as follows, $$f\left( {\mathbf{v}} \right){d^3}v = \frac{N}{\left(\pi\kappa\theta^2\right)^{3/2}}\frac{\Gamma(\kappa+1)}{\Gamma(\kappa-1/2)}\left(1+\frac{v^2}{\kappa\theta^2}\right)^{-(\kappa+1)}d^3v,$$ where $N$ is the number density, $\kappa$ is the power-law index of the kappa distribution, $\theta$ is the most probable speed, $\Gamma$ is the Gamma function.

In [2]:
%config InlineBackend.figure_format = 'retina'
import numpy as np
import matplotlib.pyplot as plt
import scipy.special as ss
import math as mt
plt.rcParams['font.size'] = 14
plt.rcParams['axes.linewidth'] = 1.5
In [3]:
v=np.linspace(-5,5,1000)
kappa=np.array([2,6,25])
plt.figure(figsize=(6,4))
for k in kappa:
    fv=k**(-3/2)*ss.gamma(k+1)/ss.gamma(k-0.5)*(1.+v**2/k)**(-k-1)
    plt.plot(v,fv,label=r'$\kappa=$%2.1f' %(k))
    
plt.plot(v,np.exp(-v**2),'k--',label='Maxwellian')
plt.yscale('log')
plt.xlim(-5,5)
plt.ylim(1e-7,1e1)
plt.legend()
plt.xlabel(r'$v/\theta$')
plt.ylabel(r'$f(v)$')
plt.show()

This plot shows an example of a Kappa distribution with several $\kappa$ parameters. If $\kappa$ is small, one can see a high energy tail. As $\kappa$ increases, the distribution approaches to a Maxwell distribution.

Generate Kappa Distribution Using Uniform Random Number

Here, I produce a Kappa distribution for particles which can be used in PIC simulation. You can check the algorithm in detail in the reference.

In [4]:
nop=100000
m=1
T=1
kappa=3
vth=np.sqrt(2*T/m)
#sig=np.sqrt(0.5*)
u1=np.random.rand(nop)
u2=np.random.rand(nop)
w1=np.random.rand(nop)
w2=np.random.rand(nop)

vx=np.sqrt(-2*np.log(u1))*np.sin(2*np.pi*w1)
vy=np.sqrt(-2*np.log(u1))*np.cos(2*np.pi*w1)
vz=np.sqrt(-2*np.log(u2))*np.sin(2*np.pi*w2)
nn=np.sqrt(-2*np.log(u2))*np.cos(2*np.pi*w2)

if(isinstance(kappa,int)):
    Ui=np.random.rand((kappa-1),nop)
    chi2=-2*np.log(np.prod(Ui,axis=0))+nn**2
else:
    Ui=np.random.rand(int(kappa-1/2),nop)
    chi2=-2*np.log(np.prod(Ui,axis=0))
    
A=np.sqrt((2*kappa-3)*T/(m*chi2))

vx=A*vx
vy=A*vy
vz=A*vz
In [12]:
ep =0.5*(vx*vx+vy*vy+vz*vz)
vv =np.sqrt(vx*vx+vy*vy+vz*vz)

plt.figure(figsize=(24,6))
plt.subplot(121)
xmin=1e-1
xmax=2e1
bins=100
x=np.logspace(np.log10(xmin),np.log10(xmax),bins)        
plt.hist(vv,bins=x,weights=1/vv,edgecolor='black',alpha=0.25)
plt.plot(x,8e3*x**2*np.exp(-x**2),'--')
#plt.plot(x,1e6*x**(-4))
plt.plot(x,2e3*x**2*ss.gamma(kappa+1)/ss.gamma(kappa-0.5)*(1.+x**2/kappa)**(-kappa-1),'k--')
#plt.legend()
plt.xscale('log')
plt.yscale('log')
plt.xlabel(r'$|v|/\theta$')
plt.ylabel('# of particles')
plt.ylim(1e-2,1e5)

plt.subplot(122)
xmin=1e-2
xmax=1e2
bins=100
x=np.logspace(np.log10(xmin),np.log10(xmax),bins)        
plt.hist(ep,bins=x,weights=1/ep,edgecolor='black',alpha=0.25)
plt.plot(x,2e4*np.sqrt(x)*np.exp(-2*x),'--')
#plt.plot(x,1e6*x**(-4))
plt.plot(x,5e3*np.sqrt(x)*ss.gamma(kappa+1)/ss.gamma(kappa-0.5)*(1.+2*x/kappa)**(-kappa-1),'k--')
#plt.legend()
plt.xscale('log')
plt.yscale('log')
plt.xlabel(r'$E/m\theta^2$')
plt.ylabel('# of particles')
plt.ylim(1e-2,1e5)
plt.show()

The left and right panels show a Kappa distribution with $\kappa=3.5$ for the speed and energy produced by using uniform random numbers. The black dashed lines are the theoretical kappa distribution function. The orange dashed lines are a fitted Maxwell distribution proposed by Oka et al. 2013. You can find that the kappa distribution generated with random numbers fits well to the theory.