Plasma via PIC (2 species – 1D)

In this page you can find a (non optimized) Matlab/Octave script devoted to solve the dynamics of charged particles in external and internal fields. The code take into account two different species (electron and ions), in order to be able to simulate large density number each particle is weighted, that mean that mass and charge are not the real one but instead are multiplied by a certain weight. Weighted particles are often called super-particles.

CODE SCHEME:

• Definition of time and space mesh
• Definition of particle number and initial condition
• Definition of external field E(x,t)
• Time cicle:
  • Evaluation of the field for each particle (nearest point)
  • Acceleration
  • Evaluation of the charge density (projected on the mesh)
  • Computation of the new field

FIELD GENERATED BY PARTICLES

Let’s start from the field generated by the charges, knowing the charge density we can solve the Poisson equation (Magnetic field due to current is neglected): \nabla^2 \bar \Phi = \bar \rho (guassian dimensionless units \bar \rho = \rho 4\pi e /(m\omega^2 )) using periodic boundary condition.

function [E,PHI] = electrostatic_1D(x,rho)

N_space = length(x);  %number of nodes
dx = x(2) - x(1);     %cell dimension

E = zeros(1,length(x)); %Electric field definition

%BC periodic, one extreme must be fixed

PHI_END = 0;

%SYSTEM MATRIX (second derivative)
e = ones(N_space,1);
A = -1/dx^2*spdiags([1*e -2*e 1*e],[-1 0 1],N_space,N_space);
A(1,end) = -1/dx^2;
A(end,end-1) = 0;
A(end,end) = -1;    %periodic condition

rho(end) = PHI_END; %periodic condition

%SOLVE POISSON EQ

PHI = A\rho';

%EVALUATION OF THE FIELD

E(2:end-1) = - (PHI(3:end) - PHI(1:end-2))/2/dx; %centered FD
E(1) = -(PHI(2) - PHI(1))/dx;
E(end) = -(PHI(end) - PHI(end-1))/dx;

end

This function is called every timestep after the definition of the overall density.

PROJECTION ON THE MESH

An important point is the projection of the quantities on the mesh, the easist procedure is the “nearest grid point approach” that mean for exemple that a particle feel the field not in its position but in the nearest grid point (we need this because the field is evaluated on a mesh but the particle can assume every position without restrictions):

function [E_sp]=field_projection(E,X_sp,x)

%FIND NEAREST NODE ID
[val, id_x] = min(abs(X_sp - x));

E_sp = E(id_x);

end

Similarly for the charge density projection.

BOUNDARY CONDITION

In this case the boundary conditions are periodic, a way to do that in Matlab/Octave is using the rem function: (Be careful, it can contain bug)

function [X_sp] = periodic_boundary(X_sp,x)

if rem(X_sp,2*x(end)) == 0 
    X_sp = sign(X_sp)*(2*x(end) - x(end)) ;
elseif rem(X_sp,x(end)) > 0
    X_sp = rem(X_sp,x(end)) - x(end);
else
    X_sp = rem(X_sp,x(end)) + x(end);
end

end

In the main an if statement check if the particle coordinate X is out from the domain if so periodic_boundary() is called. Then there are three possibilities: particle on the boundary, right to the boundary, left to the boundary.

CHARGE ACCELERATION

This section probably don’t need any explanation:

function [X_sp,V_sp] = particle(E_sp,X_sp,V_sp,dt,q,m)

%PARTICLE ACCELERATION

%EXPLICIT	
 		
 V_sp = V_sp + dt*q/m*E_sp; 
 X_sp = X_sp + dt*V_sp;	 

end