Particle Beam Example

[abigailalexjames]

Hi there,
I'm new to PIC codes and need to implement an electron beam in a plasma. I believe I need to use the injector block but I wondered if there are any workshop examples of how to create a particle beam. The documentation on the injector block just isn't clear enough and I'm not sure how to implement control of the shape and velocity of the beam.

I've seen issues/discussions of people using particle beams but not anything related to initial implementation.

I would really appreciate any help or recommendations for resources on how to do this.

Many thanks

[Status-Mirror]

Hey @abigailalexjames,

In the coming weeks I expect to be developing more resources for getting started with EPOCH, as the current documentation mainly focuses on the syntax of existing input deck blocks. This is less useful for complete beginners.

In the meantime, I'll provide a basic example here. Let me know if this is helpful or not, as it may one day become an example on the website

For our input deck, we'll create three species - two representing electrons and ions in the plasma, and a third species for injected beam electrons. Field and particle boundary conditions have been set to open, so particles can move in and out of the simulation window. However, we don't want to lose resolution for our background plasma, so we use thermal boundaries for these. Every time a macro-particle leaves the simulation window on a thermal boundary, it is replaced with an incoming particle with a temperature sampled from the initial conditions. Physically, this is supposed to represent plasma from beyond the simulation boundary pushing a macro-particle in to replace the lost one. For the beam species, load no particles initially and handle the properties in the injector block. Here I inject a 9 MeV beam with a square density profile through the x_min boundary. This simulation took me under 3 mins to run on my desktop 4-core computer, and generated 853MB of total output.

begin:control
    nx = 400
    ny = 400
    npart = 10 * nx * ny
    t_end = 100.0e-15
    x_min = 0
    x_max = 20.0e-6
    y_min = -10.0e-6
    y_max = 10.0e-6
    stdout_frequency = 10
end:control

begin:boundaries
    bc_x_min = open
    bc_x_max = open
    bc_y_min = open
    bc_y_max = open
end:boundaries

begin:species
    # Ionised hydrogen plasma
    name = Plasma_ion
    mass = 1836.2
    charge = 1.0
    frac = 0.2
    density = 1.0e24
    temp_ev = 1000.0
    bc_x_min = thermal
    bc_x_max = thermal
    bc_y_min = thermal
    bc_y_max = thermal
end:species

begin:species
    # Background plasma electrons
    name = Plasma_electron
    mass = 1.0
    charge = -1.0
    frac = 0.8
    density = 1.0e24
    temp_ev = 1000.0
    bc_x_min = thermal
    bc_x_max = thermal
    bc_y_min = thermal
    bc_y_max = thermal
end:species

begin:species
    # Electron beam
    name = Electron_beam
    mass = 1.0
    charge = -1.0
    frac = 0
end:species

begin:injector
   boundary = x_min
   species = Electron_beam
   number_density = if (abs(y) lt 3.0e-6, 1.0e10, 0.0)
   temp_ev = 1.0e3
   drift_x = 5.0e-21       # Around 9 MeV
   t_start = 0.0
   t_end = 50.0e-15
   npart_per_cell = 10
end:injector

begin:output
    dt_snapshot = 10.0e-15
    number_density = always + species
    px = always
    weight = always
end:output

I developed a simple MATLAB plotter to visualise the SDF files, and I'll attach my script below. You'll find the GetDataSDF function in epoch/SDF/Matlab, along with some other scripts required for reading SDF.

% User input
file_name = '0004.sdf';
species_name = 'Electron_beam';

% Extract data
data = GetDataSDF(file_name);
density = data.Derived.Number_Density.(species_name).data;
time = data.time;
x_edges = data.Grid.Grid.x;
y_edges = data.Grid.Grid.y;

% Create plot
x_centres = 0.5*(x_edges(2:end)+x_edges(1:end-1));
y_centres = 0.5*(y_edges(2:end)+y_edges(1:end-1));
[x_plot, y_plot] = meshgrid(x_centres, y_centres);
surf(x_plot*1.0e6, y_plot*1.0e6, density', 'EdgeColor','none');

% Figure formatting
view(2);
cbar = colorbar;
xlabel('x [\mum]');
ylabel('y [\mum]');
title(['t = ' sprintf('%2.2f',time*1.0e15) ' fs']);
ax = gca;
ax.FontSize = 16;
cbar.Label.String = 'Density [1/m³]';

which produced the following plot of the electron beam density. Here we can more clearly see the square profile in the $y$ direction for the injected beam

Let me know if you think this quick demo could be improved!

Cheers,
Stuart

[Status-Mirror]

No problem! It sounds like you want the focusing behaviour from the Gaussian beam demo, and the temporal profile from the laser pulse demo.

In EPOCH, the laser boundary theory expects lasers to be injected perpendicular to the boundary - it is always easier to rotate your target instead of having the laser at an angle. I think the only time you would need an angled injection would be if you have more than one laser in the simulation, coming in at different non-perpendicular angles.

Regarding the point in your input.deck, the laser would re-appear on the opposite side of a periodic boundary if it passed through (I had to check this recently). If the periodic boundaries are there to stop macro-particles escaping through y_min and y_max, I would instead consider open boundaries for the fields, but thermal boundaries for the particles. See the boundary deck documentation for more details.

Cheers,
Stuart

[abigailalexjames]

Hi Stuart,
Thank you. I will indeed have two laser beams which has been the main struggle so far. I have a drive beam that ideally will be at an angle followed by a probe beam perpendicular to the boundary. The initial struggle was making sure the beam is focussing correctly at an angle to accurately drive the wakefield.
It's a trade off it seems, having the drive beam at an angle makes resolving the electric field from the probe easier, however it's easier to keep the drive beam perpendicular. I originally thought I could combine the focussing gaussian beam demo with the laser at an angle demo but it turns out to be a bit trickier than that. I found a focussed beam at an angle issue here but haven't been able to combine it all so that it works just yet.

Thanks for the note about the boundaries. I felt that I should have open/thermal but it was mentioned to me that it needs to be periodic however I can't determine why. If I don't want to lose information sure however I would expect to extract the information from the domain before the lasers hit the boundaries or that the particles that leave the simulation wouldn't affect the region where I want information. (apologies for the long discussion but it's incredibly helpful to be able to ask all these questions!)

[Status-Mirror]

I've tried to get focusing Gaussian beams to come in at an angle before, but I've never quite managed it. The profiles seem to slant themselves in such a way that it doesn't look like a normal Gaussian beam if you rotate the domain back. I think you'd have more success giving a temporal profile to the rotated beam example.

[abigailalexjames]

I'll give it a try. And that is very reassuring. At least I know I'm not just making a trivial mistake!

[ajazmir]

How to load particles in the injector. I want to check x-px of Electron Beam and Plasma Electrons.
How to improve it to load particle data of electron beam from t = 0 to t=t_end.

begin:control
nx = 400
ny = 400
npart = 10 * nx * ny
t_end = 100.0e-15
x_min = 0
x_max = 20.0e-6
y_min = -10.0e-6
y_max = 10.0e-6
stdout_frequency = 10
end:control

begin:boundaries
bc_x_min = open
bc_x_max = open
bc_y_min = open
bc_y_max = open
end:boundaries

begin:species
# Ionised hydrogen plasma
name = Plasma_ion
mass = 1836.2
charge = 1.0
frac = 0.2
density = 1.0e24
temp_ev = 1000.0
bc_x_min = thermal
bc_x_max = thermal
bc_y_min = thermal
bc_y_max = thermal
end:species

begin:species
# Background plasma electrons
name = Plasma_electron
mass = 1.0
charge = -1.0
frac = 0.8
density = 1.0e24
temp_ev = 1000.0
bc_x_min = thermal
bc_x_max = thermal
bc_y_min = thermal
bc_y_max = thermal
end:species

begin:species
# Electron beam
name = Electron_beam
mass = 1.0
charge = -1.0
frac = 0
end:species

begin:injector
boundary = x_min
species = Electron_beam
number_density = if (abs(y) lt 3.0e-6, 1.0e10, 0.0)
temp_ev = 1.0e3
drift_x = 5.0e-21 # Around 9 MeV
t_start = 0.0
t_end = 50.0e-15
npart_per_cell = 10
end:injector

begin:output
dt_snapshot = 10.0e-15
number_density = always + species
px = always
weight = always
end:output