unexpected qed emission from a plasma

[ktangtar]

Hi,

My colleagues and I have been running into some weird photon emission problems. Ultimately, we ended up making the simulation really minimal to help debug. It is just a plasma with the qed module turned on. We end up seeing unexpected photons being emitted throughout the plasma.

The electron density:

The electron distribution (does not even reach mc^2):

The accumulated photon energy distribution and accumulated photon number distribution (500 fs):

A scatter plot of photons emitted with over 10keV energy:

Ex and Ey fields [V/m], the colorbar range is from the minimum to maximum cell (noise is about ten orders of magnitude lower than the sheath fields):

Is this a bug in the qed emission module? The plots are from the input deck below.

Regards, Kavin Tangtartharakul

begin:control
  nx = 100*20
  ny = 50*20

  # Size of domain
  x_min = 0 * micron
  x_max = 100 * micron
  y_min = -25 * micron
  y_max = -y_min

  field_ionisation = F

  # Final time of simulation 1000, (2640 for 800micron, 120 dt - 22 frames, change center)
  t_end = 500 * femto
  #restart_snapshot = restart_dumps0002.sdf
  stdout_frequency = 10
  use_random_seed = T
end:control

begin:qed
  use_qed = T
  use_radiation_reaction = T
  produce_photons = T
  photon_energy_min = 100*ev
  photon_dynamics = F
  qed_table_location = /scratch/08925/talia/epoch_photon_shape/epoch-4.17.16/epoch2d/src/physics_packages/TABLES
end:qed

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

begin:constant
  n_critic = 1.71e21 / cc
  z_electron = 1
  #temperature_ev = 0.3531 
  temperature_ev = 20000

  temp_SI_K = 11604 * 20000
  mass_e_kg = 9.10938e-31
  Kb_J_K = 1.380649e-23
  mKT =  mass_e_kg * Kb_J_K * temp_SI_K
  P_T_sigma = (mKT)^0.5

  dens_cm = 1.502e19 #[cm-3]
  dens_m = dens_cm / cc #[m-3], (cc = 10^-6 convert to cubic metres)
  dens_ne = dens_m
  dens_ni = dens_m / z_electron
  center_plasma = 50 * micron
  step_x = (1 / (1 + exp(-2000000 * (x - 20 * micron)))) - (1 / (1 + exp(-2000000 * (x - 80 * micron))))

  lambda0 = 0.8 * micron

  freq0 = c / lambda0
  Laser_Pd_T = 1 / freq0

  #Radius of the beam waist
  w_0 = 2.625 * micron

  #Focal length of the beam
  f_x = center_plasma

  #Rayleigh length
  Z_R = (pi * (w_0)^2) / lambda0

  #spot size evolution
  W_x = (w_0) * sqrt(1 + (f_x / Z_R)^2)

  #Radius of curvature
  R_c = f_x * (1 + (Z_R / f_x)^2)

  #Guoy phase determination
  g_phase = atan(f_x / Z_R)

end:constant

begin:laser
  boundary = x_min
  #intensity_w_cm2 = 2.9587e19 * (w_0 / W_x)
  intensity_w_cm2 = 0
  lambda = lambda0
  pol = 0.0
  phase = (2.0*pi/lambda0) * (y^2/(2.0*R_c)) - g_phase
  profile = gauss(y,0,W_x)
  t_profile = gauss(time, 36 * femto, 18 * femto)
end:laser

begin:species
  # electrons
  name = el_start
  charge = -1
  mass = 1.0 #electron mass units
  #temperature = temperature_ev
  dist_fn = (2*pi*mKT)^(-3/2) * exp(-(px^2+py^2+pz^2)/(2*mKT))
  dist_fn_px_range = (-3*P_T_sigma,3*P_T_sigma)
  dist_fn_py_range = (-3*P_T_sigma,3*P_T_sigma)
  dist_fn_pz_range = (-3*P_T_sigma,3*P_T_sigma)
  number_density = dens_ne * step_x
  number_density_min = n_critic * 1.0e-6
  npart_per_cell = 4
  identify:electron
end:species 

begin:species
  name = photons
  charge = 0.0
  mass = 0.0
  npart = 0.0
  number_density = 0.0
  identify:photon
end:species

begin:species
  # ions
  name = ions
  charge = 1
  mass = 1835.16 #electron mass units, Hydrogen
  #mass = 359420 #electron mass units, Gold
  temperature = 0
  #density_min = n_critic /(1000* z_electron)
  number_density = dens_ni * step_x
  number_density_min = n_critic * 1.0e-6
  npart_per_cell = 4
  immobile = T
end:species

begin:subset
  name:high_energy_el
  include_species:el_start
  gamma_min = 2
end:species

begin:subset
  name:all_el
  include_species:el_start
end:species

begin:subset
  name:photons_sub
  include_species:photons
end:species

begin:dist_fn
  name=xy
  ndims=2
  dumpmask=always
  direction1=dir_x
  direction2=dir_y
  include_species:el_start
end:dist_fn

begin:output
  name=dist
  file_prefix = dist
  dt_snapshot = 100 * femto
  distribution_functions = always + single + snapshot
end:output

begin:output
  name = charge_dens
  file_prefix = charge_dens

  # Number of timesteps between output dumps
  dt_snapshot = 100 * femto
  dt_average = 4 * Laser_Pd_T

  charge_density = always + single

end:output

begin:output 
  name = absorp_totalenergy
  file_prefix = fields_var

  # Number of timesteps between output dumps
  dt_snapshot = 100 * femto
  dt_average = 4 * Laser_Pd_T

  absorption = always
  total_energy_sum = always + species

end:output

begin:output 
  name = fields_var
  file_prefix = fields_var

  # Number of timesteps between output dumps
  dt_snapshot = 100 * femto
  dt_average = 4 * Laser_Pd_T

  # Properties on grid
  grid = always
  ey = always + average + snapshot + single
  ex = always + average + snapshot + single
  bz = always + average + snapshot + single

end:output

begin:output

  name = regular_output_particles
  file_prefix = particles_properties

  # Number of timesteps between output dumps
  dt_snapshot = 25 * femto

  # Properties at particle positions

  particles = high_energy_el + single
  particle_energy = high_energy_el + single
  work_x_total = high_energy_el + single
  work_y_total = high_energy_el + single
  particle_weight = high_energy_el + single
  px = high_energy_el + single
  py = high_energy_el + single
end:output

begin:output
  name = all_regular_output_particles
  file_prefix = all_particles_properties

  # Number of timesteps between output dumps
  dt_snapshot = 25 * femto

  # Properties at particle positions

  particles = all_el + single
  particle_energy = all_el + single
  work_x_total = all_el + single
  work_y_total = all_el + single
  particle_weight = all_el + single
  px = all_el + single
  py = all_el + single
end:output

begin:output
  name = restart_dumps
  file_prefix = restart_dumps
  dt_snapshot = 500 * femto
  restartable = T
  dump_first = F 
end:output

begin:output 
  name = current_var
  file_prefix = current_var

  # Number of timesteps between output dumps
  dt_snapshot = 100 * femto

  # Properties on grid
  grid = always
  jx = always + single
  jy = always + single

end:output

begin:output
  name = regular_output_photon
  file_prefix = photon_properties

  # Number of timesteps between output dumps
  dt_snapshot = 25 * femto

  # Properties at particle positions

  particles = photons_sub + single
  particle_energy = photons_sub + single
  work_x_total = photons_sub + single
  work_y_total = photons_sub + single
  particle_weight = photons_sub + single
  px = photons_sub + single
  py = photons_sub + single

end:output

[Status-Mirror]

Hi @ktangtar,

I don't think this is self-heating - from small-scale test simulations, I think your electron temperature stays stable over the simulation run-time. I'm pretty sure you're not gaining and radiating non-physical energy. I've included my tiny parameter test deck at the end of this message.

In my tiny test, the electrons reached a peak energy of around 100 keV. It's worth pointing out that all electrons have the ability to radiate when you switch on QED - not just the ones with KE over the rest mass energy 511 keV. All electrons would then start recoiling from sampled photon emissions, and if those photons have energy over the value set by your photon_energy_min key, then they are also added to the grid as macro-photons.

The maximum possible photon energy is the electron KE, and your photon_energy_min is 100 eV. Hence, I would expect your macro-photon distribution to span from 100 eV to ~100 keV, which is roughly what I see here. Even though your fields inside the plasma are small, they are still non-zero and capable of producing a small amount of radiation. If there was a strong laser in the simulation, I think you'd find that this plasma emission was negligible in comparison.

If that emission is problematic, you can reduce the field noise by increasing the particles per cell, compiling the code with a higher order particle shape, or running the code with smooth_currents = T in the control block.

Hope this helps, Stuart

begin:control
  nx = 10
  ny = 10
  t_end = 500.0e-15
  x_min = 0
  x_max = nx * 50.0e-9
  y_min = 0
  y_max = ny * 50.0e-9
  stdout_frequency = 10
end:control

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

begin:species
  name = proton
  npart_per_cell = 4
  density = 1.5e25
  temp_ev = 0.0
  immobile = T
  identify:proton
end:species

begin:species
  name = electron
  npart_per_cell = 4
  density = 1.5e25
  temp_ev = 20.0e3
  identify:electron
end:species

begin:output
  dt_snapshot = t_end / 10.0
  grid = always
  px = always + species
  py = always + species
  pz = always + species
  weight = always + species
  temperature =  always + species
end:output

[ktangtar]

Thank you for the reply @Status-Mirror.

I think for these slightly relativistic emitting electrons (1 mc^2), in the weak noise fields (1e10 V/m), it should be extremely unlikely for the photons to have energies anywhere near that of the electrons. The critical photon emission energy should be less than 1 eV?

Using the input deck from the original post, and turning on the laser (I=3e19 W/cm^2, E0=1.5e13 V/m), the photon's energy and number distributions look like:

Like you mentioned in your comment, we should expect that the emission with the laser will be much more frequent and energetic than the background noise emission of the photons. The fields of the laser are ~1000 times that of the noise and we bring electrons to energies 10 times higher.

[Status-Mirror]

Few quick things to consider: I wouldn't think of 3e19 W/cm2 as a very high intensity for qed emission. This kind of radiation is parametrised by $\eta^2$,

$\eta^2 = \frac{1}{E_s^2} \left( \left| \gamma \pmb{E} + \frac{1}{m_e c} (\pmb{p} \times c\pmb{B})\right|^2 - \frac{1}{(m_e c)^2}(\pmb{E}\cdot\pmb{p})^2 \right)$

and EPOCH uses a look-up table to get emission probabilities and photon energies. I believe these tables only have values down to $\eta^2=10^{-5}$, and any $\eta^2$ less than this is assumed to be $10^{-5}$. This is what the truncation warning is referring to when you run the code.

Maybe this is why you get radiation without the laser, and why the distribution is similar for 3e19 - this probably won't give you a high enough $\eta^2$ to climb far up the lookup table. Normally this is fine, as you get significantly higher emissions from the higher intensities.

To get a strong qed emission and noticeably higher $\eta^2$, I would say you'd need closer to 1e21W/cm2.

It may help to estimate this $\eta^2$ parameter for your system, to see if this is actually the cause of the weird behaviour.

Hope this helps, Stuart

[ktangtar]

@Status-Mirror Yeah, I think you bring up a good point about the truncated tables used in the probabilistic emission algorithm. With the laser, eta^2 probably mostly lower than 10^-5.

If this is the main source of the problem, then I guess it becomes pretty problematic when you are simulating something emitting at the lower end of the tables. Because all the particles with minuscule emission parameters get pushed up significantly to the lowest levels of the tables? Perhaps it makes more sense to have no emission for those or make the tables larger?

[Status-Mirror]

Hey @ktangtar

CORRECTION: I believe $\eta$ goes down to $10^{-5}$, not $\eta^2$ as I previously said.

On table extension, we would run into a problem. The QED model is based on the work of Ridgers et al (2014), which asserts the probability of emission is a function of $\eta$ only. This relies on the quasistatic assumption, which forces $a_0 \gg 1$, where $a_0$ is the laser pump strength (normalised vector potential):

$a_0 = \frac{eE_0\lambda}{2\pi me c^2} \approx 8.5\times 10^{10} \sqrt{I{W/cm^2} \lambda_{\mu m}^2}$

For a 1 $\mu m$ laser at $3\times 10^{19} \text{ Wcm}^{-2}$, you have $a_0 \approx 5$, which may introduce approximation issues. The model would break down at lower values of $\eta$, and different quantities would need to be considered - this would require a re-write of the radiation module.

Perhaps in your case, it would be more sensible to cut off radiation at $\eta=10^{-5}$. To achieve this quickly, locate the file hsokolov.table stored in _epoch/epoch2d/src/physicspackages/TABLES

The second line currently reads:

-5 0.7189960759048358

To disable emission for $\eta < 10^{-5}$, modify this line to:

-5 1.0

In case you're curious, the two columns in this table describe $\log{10}\eta$, and $\log{10}h(\eta)$, where $h(\eta)$ is defined in Appendix A of Ridgers et al (2014). Changing the log term to 1 switches $h(\eta)$ to zero, and also sets the photon emission rate to zero according to (2) in Ridgers et al (2014).

There are also copies of this table for epoch1d and epoch3d, so change whichever is relevant. If this helps, we may add a new input deck key to allow users to switch off emissions for low $\eta$ without manually editing the QED table files. We also anticipate the QED routines to be updated by the University of York group in the coming year or two.

I hope this was helpful, Stuart

[ktangtar]

@Status-Mirror Wow, thanks, I think this will be really helpful and we will give it a try! Will let you know our results.

[ktangtar]

@Status-Mirror The emission rate is proportional to $\eta h(\eta) $. Since we want particles with really low $\eta$ to not emit, -5 1.0 ($\eta=10^{-5},~h(\eta)=10$) does not work. Instead, we tried -5 -1e16 ($\eta=10^{-5},~h(\eta)=0$). We have tested this for the simulation with just a plasma, and it is able to remove all the "noise" emission we were seeing before. Do you think there may be some unintended consequences in other parts of the code for doing -5 -1e16?

[Status-Mirror]

Yes you're right, I got confused with my numbers and mixed them up. The look-up tables give $\log{10}\left(h(\eta)\right)$, and I told you to set this to 1, as if you were setting $h(\eta)$ to 1 to make $\log{10}\left(h(\eta)\right)=0$, when you actually want to set $\log_{10}\left(h(\eta)\right)\rightarrow - \infty$ to make $h(\eta)\rightarrow 0$ as you do with your -1e16.

On unintended consequences, I would be surprised if there were any. The hsokolov lookup table is only used in the delta_optical_depth subroutine, which updates the travelled "optical depth" of an electron, $\tau$. Once $\tau$ reaches a pre-calculated, randomly-sampled depth of emission, an emission event is calculated and $\tau$ is reset. The only effect this should have on the code is to stop low $\eta$ contributions to $\tau$, removing the low field emissions as you see. Photon energies are sampled from a different table, even if emission in a low $\eta$ field occurred. For electron trident photon emission and Breit-Wheeler pair production, (both in pairprod.table), it looks like this had already been done by setting their equivalent table entries to -50.

The only point to remember is that we may have switched off erroneous emissions, but we still aren't calculating the correct emission for $\eta$ in this region. I would just assume this is small compared to emission in laser fields.

[Status-Mirror]

The QED modelling group from the University of York has submitted an EPOCH extension which switches to an analytic expression for low $\eta$. This has now been merged into EPOCH, and should hopefully resolve this problem.