Warwick-Plasma / epoch

Particle-in-cell code for plasma physics simulations
https://epochpic.github.io
GNU General Public License v3.0
185 stars 59 forks source link

Laser group velocity is greater than light speed with field_order=4 #647

Closed XiangyanAn closed 7 months ago

XiangyanAn commented 8 months ago

Hello, I just found that with field_order=4 in the control block, the group velocity of a laser pulse in the simulation could be greater than the speed of light. I'm using the release epoch-4.19.0, which is unzipped as epoch-4.18.0-2022-12-14, and taking 1d simulation. The input.deck is :

begin:constant lambda0 = 800e-9 k0 = 2 pi / lambda0 omega0 = c k0 tau0 = 2 * pi / omega0

a0 = 1.4 tau_laser = 8.0 lambda0 / c delay = tau_laser 3.0

end:constant

begin:control x_min = -30.0 lambda0 x_max = 30.0 lambda0

nx = 60 * 30

t_end = delay + (0 - x_min) / c + 4000.0 * tau0

field_order = 4 end:control

begin:boundaries bc_x_min = simple_laser bc_x_max = simple_outflow end:boundaries

begin:window move_window = T window_v_x = c window_start_time = (x_max - x_min) / c bc_x_min_after_move = simple_outflow bc_x_max_after_move = simple_outflow end:window

begin:laser boundary = x_min amp = a0 me omega0 * c / qe lambda = lambda0 polarisation_angle = 0.0 phase = 0.0 t_profile = gauss(time, delay, tau_laser) end:laser

begin:output dt_snapshot = t_end / 50 full_dump_every = 0

grid = full ex = full ey = full ez = full bx = full by = full bz = full end:output

After simulation, the laser pulse electric fields of the files 0001.sdf and 0050.sdf are, respectively,

ey-0001

ey-0050

Since the window is moving at the light speed, we can see that the pulse is faster than light in vacuum. Moreover, I also calculated the group velocity according to the centroid position of the laser pulse: $x_c = \int x E_y^2\mathrm d x/ \int Ey^2\mathrm d x$. The difference of the centroid position and time of 0001.sdf and 0050.sdf also gives, $(x{c,50}-x{c,1})/[c(t{50}-t_{1})]=1.0036>1$.

I'm not sure whether it is a feature of the fourth field order, or a bug. Thans for any help.

Status-Mirror commented 7 months ago

Hey @XiangyanAn,

Are you sure there's a problem here? When I calculate your t_end, I get around 10.8ps. In 10.8ps, the start of the laser can travel around 3.24mm - which is beyond your window. This is what we observe, isn't it?

Cheers, Stuart

XiangyanAn commented 7 months ago

Hi @Status-Mirror ,

We can notice that there is a delay of $t_d=24\tau0$ for the laser pulse peak and it is injected at the $x\text{min}$ boundary. Therefore, the pulse peak will arrive at the position $x=0$ at the time $t=td+\left|x\text{min}\right|/c$. Then the laser pulse will travel for $4000\tau_0=10.674 \text{ps}$ and it will arrive at 3.2 mm at the simulation end, which should be still in the window.

However, as we observed, it is not. I also calculated the group velocity quantitatively in the previous comment, which is also greater than the light speed.

Thanks, Xiangyan An

Status-Mirror commented 7 months ago

Hey Xiangyan An,

Sorry, I misread your initial input deck - I was assuming your x_min was 0. I can confirm you do have a laser which has travelled too far. I was playing around with some values, and I think the issue is due to resolution. When I double the cell-count, I find the laser travels less distance. The second-order field solver seems less sensitive to the cell-size, so this may be better for vacuum simulations.

resolution_issues

The 4th-order field solver calculates the gradients for Maxwell's equation using more cells. In a noisy plasma environment, this can smooth out the spatial field gradients and reduce noice. It seems this method is doing more harm than good in your vacuum example though.

Hope this helps, Stuart

TomGoffrey commented 7 months ago

Numerical dispersion can result in group (and phase) velocities greater than c (or indeed less than). I would suggest reading "A systematic approach to numerical dispersion in Maxwell solvers" by Blinne et al (2018) as a starting point to understand more.

XiangyanAn commented 7 months ago

Hi @TomGoffrey,

Thanks for your help. I used to think that numerical dispersion would only result in group velocities less than $c$. Otherwise, there might be causality problems. I have also observed group velocity less than $c$ with field_order = 2. That's why I'm surprised by the group velocity greater than $c$ with field_order = 4. I now have a deeper understanding of numerical dispersion. Thanks again for your reminder.

Also thank @Status-Mirror for the further detailed check. This problem can indeed be alleviated by using more cells.

amarsarari commented 7 months ago

Hey Xiangyan An,

Sorry, I misread your initial input deck - I was assuming your x_min was 0. I can confirm you do have a laser which has travelled too far. I was playing around with some values, and I think the issue is due to resolution. When I double the cell-count, I find the laser travels less distance. The second-order field solver seems less sensitive to the cell-size, so this may be better for vacuum simulations.

resolution_issues

The 4th-order field solver calculates the gradients for Maxwell's equation using more cells. In a noisy plasma environment, this can smooth out the spatial field gradients and reduce noice. It seems this method is doing more harm than good in your vacuum example though.

Hope this helps, Stuart

amarsarari commented 7 months ago

Hi sir can you please provide me matlab code to see this electric field Ey vs x

Status-Mirror commented 7 months ago

Hey @amarsarari,

First I open MATLAB, find the EPOCH directory in the internal file explorer, and right click to select "Add to path -> Selected folders and subfolders". This adds all the SDF files from the EPOCH directory to the path, so we can open them.

Once I've done this, I run this script:

% plot_Ey.m
% Plots the Ey field from a specific SDF file of a 1D EPOCH simulation

%% User input

file_name = "0050.sdf";

%% Extract data

% Pull variables from SDF file
data = GetDataSDF(file_name);
Ey = data.Electric_Field.Ey.data;
x_edges = data.Grid.Grid.x;

% Format variables
x_centres = 0.5*(x_edges(2:end) + x_edges(1:end-1));

%% Create plot

plot(x_centres * 1.0e6, Ey, 'LineWidth', 2);
box on;
grid on;
xlabel('x [\mum]');
ylabel('Ey [V/m]');
ax = gca;
ax.FontSize = 16;

For the multiple lines shown in my figure, you can repeat this process and stack the lines by running hold on, and re-running this script with a new filename. A legend can be added using the legend("Line 1 name", "Line 2 name", "Line 3 name") command.

Hope this helps, Stuart