Betatron radiation

[StevE-Ong]

Hi,

I am simulating a betatron radiation from LWFA electron beam. The radiation param I used is as attached. I run a test for 2000 steps (from 38000 to 40000) and plot the spectrum. The result I obtained is kind of weird, especially the angular part. The script I used for plotting is from https://picongpu.readthedocs.io/en/0.5.0/usage/plugins/radiation.html#pausch2014. The N_observer = 256 but the resulting angular part is divided into larger chunk as in the second figure. I am running this on 32 K40 GPUs on ARIS and it took 13hours. I hardly find a reference using picongpu on betatron radiation and most of them are on Thomson scattering. Is there any advises on this? Thanks.

radiation_param.zip

On the other hand, is there any method to plot the angular distribution of the radiation, i.e. distribution over theta and phi?

[sbastrakov]

@PrometheusPi @Anton-Le are most familiar with the radiation plugin. Could you take a look?

[PrometheusPi]

Please excuse the late reply @StevE-Ong. You are right, we mainly use the radiation plugin for (nonlinear) Thompsons scattering since computing radiation with the Liénard–Wiechert potential approach we use in the radiation plugin is quite expensive and for betatron radiation a simpler synchrotron radiation model is quite suitable (but of course less correct).

The frequency definition you chose appears to be valid - you go from 0.01keV to 10keV on a logarithmic scale.

For the observation directions, you chose a pattern that is similar to a quarter dome observation direction we use for the Kelvin Helmholtz example. A 16th dome (0.25pi solid angle) is quite a large field of view for forward-directed radiation. Your angle go from phi in [-pi/8, +pi/8] and theta in [-pi/8, +pi/8] with phi iterating through all of its 16 values before increasing theta. Your observation direction is defined as

vector_64( sinTheta*cosPhi , cosTheta, sinTheta*sinPhi )

Thus at index 0 you are looking into phi=-pi/8 and theta=-pi/8 into (-0.15, +0,92, -0.36), at index 128 you are looking into phi=-pi/8 and theta=0.02 and so on... Thus the chunks you see are the fixed theta regions. They are chunks with a close to fixed maximum frequency due to the Nyquist limit we apply on all spectra to ensure to only make predictions on frequencies we were able to resolve in retarded time. So far, we either fixed either theta of phi to make a 2D plot as above or we used paraview to put frequency and observation direction onto a hollow sphere. However, for the later there is no fixed workflow. If you want to visualize it such a way, I would recommend defining a 3D numpy array with a high enough resolution to translate the radiation data into a sphere.

Regarding your data, the "boxes" from fixed theta seem to overlap for y-axis > 0. This seems to be not correct. Did you alter observation directions during the simulation by recompiling picongpu?

[StevE-Ong]

Thanks @PrometheusPi for your replies. I don't quite understand the angle part. The rad_data_viaModule.get_Spectra().shape has the shape (256, 2048) only, which means that there is only one angle theta is observable while phi is fixed? or does it sum over phi for different theta? To find the relation of the index and angle, I made the following loop:

N_observer = 256
N_angle_split = 16

N_theta = N_observer / N_angle_split

angle_theta_start = -np.pi/8
angle_theta_end = np.pi/8

angle_phi_start = -np.pi/8
angle_phi_end = np.pi/8

for i in range(0,256,1):
    for j in range(0,1,1):
        delta_angle_theta =  (angle_theta_end - angle_theta_start) / (N_theta-1.0)
        delta_angle_phi =  (angle_phi_end - angle_phi_start) / (N_angle_split-1.0)
        theta = i * delta_angle_theta + angle_theta_start
        phi = j * delta_angle_phi - angle_phi_start
        sinTheta = np.sin(theta)
        cosTheta = np.cos(theta)
        sinPhi = np.sin(phi)
        cosPhi = np.cos(phi)
        print(i,j,theta,phi)

it gives

0 0 -0.39269908169872414 0.39269908169872414
1 0 -0.34033920413889424 0.39269908169872414
....
8 0 0.026179938779914924 0.39269908169872414
...
128 0 6.309365245959501 0.39269908169872414

If this is correct, then index 8 would be nearest to the axis. It seems that the current set up doesn't have a spectrum fixed at phi=0?

Regarding the "boxes", I did not do anything during the run. I have been getting this constantly, even I change from log to linear frequency. This is another example for using linear frequency I plot as follows:

from picongpu.plugins.data import RadiationData
from mpl_toolkits.axes_grid1.inset_locator import inset_axes

# access HDF5 radiation file
radData = RadiationData(path_to_sim+"simOutput/radiationHDF5/e_radAmplitudes_{}_0_0_0.h5".format(ts))

# get frequencies
omega = radData.get_omega()

# get all observation vectors and convert to angle

vec_n = radData.get_vector_n()

theta = np.arctan(vec_n[:, 0]/vec_n[:, 1]) 
# get spectrum over observation angle
spectrum = radData.get_Spectra()

# plot radiation spectrum
fig, ax = plt.subplots()

img = ax.pcolormesh(
    omega/omega_laser, 
    theta*1000, 
    spectrum, 
    norm=LogNorm()
)
ax.set(ylabel=r"$\theta \, \mathrm{(mrad)}$",
       xlabel=r"$\omega_\mathrm{rad}\,(\omega_\mathrm{L})$",
      )
cbar = plt.colorbar(img)
cbar.set_label(r"$\frac{\mathrm{d}^2 I}{\mathrm{d} \omega \mathrm{d} \Omega} \, \mathrm{(Js)}$")
plt.tight_layout()
plt.show()

The corresponding spectrum at index 8 is where omega_coh is the boundary between incoherent and coherent radiation. This is somehow close to the spectrum in FIG. 3 in the paper Nuclear Inst. and Methods in Physics Research, A 909 (2018) 419. The expression for omega_coh in this paper, omega_coh=c/(2 ne^{1/3}) doesn't seem to have the correct unit for frequency. Could you comment on this? The plasma density I use here is ne=4.3e18 cm^-3 with laser wavelength 800nm.

[PrometheusPi]

@StevE-Ong Neither is correct. The axis with 256 entries is the index of observation. In your case, since you defined two angles in radiationObserver.param, the index is just a mapping [0, 256-1] to phi and theta values which describe observation angles. To plot a line scan, you could e.g. fix a specific theta or phi value and plot the 16 values that are part of the remaining subset.

Your python code contains a slight mistake. It should read:

N_observer = 256
N_angle_split = 16

N_theta = int(N_observer / N_angle_split)

angle_theta_start = -np.pi/8
angle_theta_end = np.pi/8

angle_phi_start = -np.pi/8
angle_phi_end = np.pi/8

delta_angle_theta =  (angle_theta_end - angle_theta_start) / (N_theta-1.0)
delta_angle_phi =  (angle_phi_end - angle_phi_start) / (N_angle_split-1.0)

for i in range(0,N_theta,1):
    for j in range(0,N_angle_split,1):
        theta = i * delta_angle_theta + angle_theta_start
        phi = j * delta_angle_phi - angle_phi_start
        sinTheta = np.sin(theta)
        cosTheta = np.cos(theta)
        sinPhi = np.sin(phi)
        cosPhi = np.cos(phi)
        print(i,j,theta,phi)

# alternative with a single index as done in PIConGPU
for index in range(0,N_observer,1):
    i = index // N_angle_split
    j = index % N_angle_split
    theta = i * delta_angle_theta + angle_theta_start
    phi = j * delta_angle_phi - angle_phi_start
    sinTheta = np.sin(theta)
    cosTheta = np.cos(theta)
    sinPhi = np.sin(phi)
    cosPhi = np.cos(phi)
    print(i,j,theta,phi) 

This the angles are not correct in your output. The closest to the y-axis is at i = 7 or i = 8 (both have the same distance to the axis). Or at index = 112 till index = 143

The error with the boxes might originate from the mapping to a single angle:

vec_n = radData.get_vector_n()
theta = np.arctan(vec_n[:, 0]/vec_n[:, 1]) 

@StevE-Ong Could you provide a list based on

for i in range(256):
print(i, " : ", vec_n[i, :], theta[i])

Perhaps that helps me understand your mapping.

[StevE-Ong]

@PrometheusPi The list for angle mapping is here: angle_list.txt

[PrometheusPi]

@StevE-Ong Thanks a lot - I will have a look.

[PrometheusPi]

Okay, the mapping of a 3D vector that covers a 2D position on the unit sphere onto a 1D angle causes the strange overlapping structures.

To illustrate that, I plotted all consecutive theta values according to their index:

Thus if you plot the following image

test_rad = np.zeros((256, 512))
for i in range(256):
    test_rad[i, :] = i/256. # new color for every horizontal line
    test_rad[i, 512-np.abs(i-128):] = np.NaN # similar cutoff as in radiation plugin via Nyquist limit

plt.pcolormesh(np.arange(512), np.arange(256), test_rad)

you get this simple structure:

If I now use the theta values from your list, I get the following image: The asymmetric pattern come from the overlapping color regions, since pcolormesh seems to overlap regions previously drawn if indices jump back.

To get a single line covering the forward direction I would suggest plotting:

 plot_angle = np.linspace(-np.pi/8., +np.pi/8., 16)
 img = ax.pcolormesh(
    omega/omega_laser, 
    plot_angle, 
    spectrum[128:144, :], 
    norm=LogNorm()
    )

I hope this input helps you.

[StevE-Ong]

Thank you @PrometheusPi

The overlapping doesn't appear over the entire range now. The mesh boxes along omega may be due to the resolution, which is not the problem for now.

Having some axis limit and log scale, it looks like this

which agrees with 1D spectrum at index 128

Regarding the omega_coh in the publication, is the expression correct? I am using the expression omega_coh = c * ne^{1/3} from another publication (FIG 1 in https://journals.aps.org/pre/abstract/10.1103/PhysRevE.92.023305). I am not sure how the expression come from, but I am appreciate if you could let me know any references for that. Thank you.

[PrometheusPi]

You are welcome @StevE-Ong.

The boxes in omega are common for limited resolution. You are totally right, they are no problem. You plots look great and correct to me.

Regarding omega_coh - this appears to be correct to me. I assume that you use a form factor other than radFormFactor_incoherent or radFormFactor_coherent, thus the drop from 1e-22 to 1e-26 you see is the transition from coherent to incoherent radiation as described in the paper by A. Gonoskov et al. you referenced. Please note, Fig. 1 in that paper focuses on photon energy (hbar * omega) thus it multiplies the spectral intensity with omega which is not done in your plot, thus you see a slowly decreasing incoherent region at frequencies above omega_coh.

(Please be aware that the radiation plugin in PIConGPU is classical and will deviate from the predictions in A. Gonoskov et al. in cases of chi>1 and also will have no energetic cutoff.)

[PrometheusPi]

If you want radiation prediction based on the models in A. Gonoskov et al. please use the synchrotron model (which is currently only extremely sparsely documented: https://picongpu.readthedocs.io/en/0.5.0/models/photons.html)

[StevE-Ong]

Multiplying omega look very close to FIG. 1 in A. Gonoskov et al. So far the radiation calculate here is betatron radiation and mostly classical. Since you mentioned the synchrotron model, is it valid also for betatron type radiation? I often see is the calculation for photon generation or QED emission, but I am skeptical about its validity for betatron radiation. But I could give it a try and I believe it would be faster than the radiation plugin.

[PrometheusPi]

It is definitely faster than the radiation plugin (but it might also run into out-of-memory problems due to the high number of photons). I would it consider valid for betatron radiation if you are not looking for spectral peaks as in Thompson scattering.

[StevE-Ong]

Alright...I will give it a try and come back for an update. Thank you @PrometheusPi for your time.

[StevE-Ong]

Hi,

I have recalculated the radiation spectrum of the same problem by changing the range of energy and observation angle. We make some analytical prediction to omega_coh, omega_critical and the peak of synchrotron radiation. The critical frequency is around 2.9 keV. Then I plot the spectrum as follows

The position of the lines seems reasonable. The first line separates between coherent and incoherent radiation. The second line is the peak of synchrotron spectrum, which is 0.29*critical frequency (the third line). What confuses me is the harmonic at the higher frequency. I did see some harmonic from previous publication of picongpu for radiation from LWFA. I expect to see some harmonic, and they end at critical frequency. I am not sure what is this. Is there something not correct with my params or the way I analyze has problems? Thanks.

radiation_param.zip

[PrometheusPi]

Hi @StevE-Ong, we observed these harmonics during some first betatron radiation tests too. I doubt that they are physics-wise correct and assume that the originate from floating point accuracy or discretization. (You could cross-check this by switch precisions from 32bit to 64 bit, but the radiation kernel performs internal casts that might hide these effects even after switching.)

The harmonics you might have seen came from nonlinear Thomson scattering. As long as you have no (high amount of) relativistic electrons with a similar gamma of ~30 passing head-on through the laser, these should not be harmonics from Thomson scattering.

[StevE-Ong]

Hi @PrometheusPi, you mean the precision in precision.param? It seems the precision is for the special operations. So by default picongpu compute everything with 32bit?

[PrometheusPi]

Yes, I meant switch from 32bit to 64bit in precision.param. Yes, by default everything except some numerically critical parts are 32bit.

[StevE-Ong]

Hi, I resumed the calculation with 64bit in precision.param. Indeed, the harmonics disappeared. The shape of the incoherent radiation also changed.

PARAM_PRECISION precision64Bit

PARAM_PRECISION precision32Bit

[StevE-Ong]

I have a question regarding the coherent part in the spectrum. In this paper ,

The scattering of laser light is coherent while the betatron radiation in the x-ray frequency range [6] is incoherent.

In this simulation, we did not see any laser scattering (particle energy filter was used). So I suppose this coherent radiation comes from the radiation for those particles where their inter-particle distance is less than the radiation wavelength, or it comes from the particles inside the macro-particle. Or what is the correct understanding?

[sbastrakov]

@PrometheusPi tagging just to make sure

[PrometheusPi]

Thanks @sbastrakov for tagging. The light scattered from the laser only scatters from particles in the laser field. These particles (in LWFA) do not reach energies much above the normalized laser field strength a_0: beta * gamma ~ a_0. Thus if your energy (gamma) filter is set as gamma >> a_0, particles in the laser field are ignored. This is the common case for today's laser field strength in LWFA and gamma filters of 10.

Might this explain, why you do not see laser scattered light in your setup with particle filters @StevE-Ong?

[StevE-Ong]

Thanks @PrometheusPi. I think this might explain that the coherent part (below \omega_coh) is the radiation from the electrons in the macroparticle?

[PrometheusPi]

@StevE-Ong Most likely yes. You can access the shape function for the grid you used by comparing the coherent enhancement you see with the form-factors equation defined here.

[PrometheusPi]

@StevE-Ong We fixed a bug in the radiation plugin that caused additional spectral modulations at high frequencies #4231. Perhaps you want to rerun your simulation and check whether the high frequency spectrum changes.