Skip to content

Tutorial 2b: Generating the Synchrotron Matrices


With the Compton matrices generated, we can move on to the Synchrotron matrices.

Generating the Emission Matrices

Particle Names

Synchrotron radiation is an emissive process 123 that unlike binary processes includes only three particles: the incoming particle (electron) that will undergo some emissive process in the presence of some external field (magnetic field) to after which it will remain along with a produced particle (photon). Therefore for the emissive synchrotron interaction of electrons can declared as:

julia
using Diplodocus

name1 = "Ele"
name2 = "Ele"
name3 = "Pho"
type = "Sync"

where type indicates the type of emissive interaction (see Implemented Collisions).

Momentum Space Grids

The momentum grids are just the same as the Compton matrix setup:

julia
    p_low_Ele = -3.0
    p_up_Ele = 9.0
    p_grid_Ele = "l"
    p_num_Ele = 288
    u_grid_Ele = "u"
    u_num_Ele = 1
    h_grid_Ele = "u"
    h_num_Ele = 1

    p_low_Pho = -16.0
    p_up_Pho = 9.0
    p_grid_Pho = "l"
    p_num_Pho = 600
    u_grid_Pho = "u"
    u_num_Pho = 1
    h_grid_Pho = "u"
    h_num_Pho = 1

Define External Parameters

For synchrotron radiation we must tell the code what magnetic field strengths to consider. This is done by defining the Ext vector, which for type=Sync is taken to be a vector of magnetic field strengths in SI units, i.e. Tesla. Here we will again be generous and ask the code to sample a wider range of magnetic field strengths than we expect to use in this tutorial so that we can re-use these emission matrices for future simulations. The range we will sample is from 1e-10\,\mathrm{T}=10^{-6}\,\mathrm{G} to 1e1\,\mathrm{T}=10^{5}\,\mathrm{G}:

julia
    Ext::Vector{Float64} = [1e-10,1e-9,1e-8,1e-7,1e-6,1e-5,1e-4,1e-3,1e-2,1e-1,1e0,1e1];
nothing #hide

Monte-Carlo Setup

For emissive processes there are fewer incoming particles involved, so we can sample far more incoming points:

julia
    numLoss = 256
    numGain = 256
    numThreads = Threads.nthreads()

WARNING

By default this script uses the maximum number of available threads Threads.nthreads(), if you would prefer the scripts to run with fewer threads then modify the variable numThreads to the desired amount.

Unlike binary collisions, we do not need to deploy the importance sample to improve the MC accuracy of emissive processes, this is because there are not dependent variables in the outgoing states, so all states get sampled equally. But for future-proofing, we can still define the scale factors to sample. Here we will just make it sample scale=0.0:

julia
scale = 0.0:0.1:0.0

Where to Save the mission Matrices

julia
fileLocation = joinpath(pwd(),"Data")

Final Setup and Running

We can bundle all the above user defined parameters into a tuple called Setup. This is generated by the UserEmissiveParameters() function, which reads the above parameters from Main. It also generates the fileName under which the emissive matrices will be saved and stored.

julia
(Setup,fileName) = UserEmissiveParameters()

WARNING

It is not advised to edit fileName as its format is used internally by DiplodocusCollisions and DiplodocusTransport to identify the appropriate collision matrices to load.

Finally, we can run the Monte-Carlo integration using the function EmissiveInteractionIntegration:

julia
EmissiveInteractionIntegration(Setup)

Checking the Collision Matrices

Unlike binary processes, particle number is not conserved in an emissive process. Energy should be conserved but currently the change of the electron distribution for synchrotron radiation is included as an analytic force term (see the tutorial 2c.) as such we cannot examine the energy conservation from the emission matrices alone.

Plotting the Spectra

What we can do though is visualise the resulting photon spectrum using the plotting function InteractiveEmissiveGainLossPlot from DiplodocusPlots. This creates an interactive plot in which the user can look at different incoming states and view the resulting outgoing spectra. To use this function, first load the file containing the collision matrices using EmissionFileLoad_Matrix:

julia
Output = EmissionFileLoad_Matrix(fileLocation,fileName);
InteractiveEmissiveGainLossPlot(Output)

This will display the plot in a separate window where you can interact with it. There are sliders for the the discrete bins of the incoming particle's (electron's) momentum p1, and their propagation directions, in spherical polars u (u=cosθ) and h (phi). There are then also sliders for the outgoing particle's (photon's) direction with the spectra of outgoing states as a function of momentum being plotted.

Now on to the simulation itself!