n-Photon Compton Scattering Process

In this file, we set up an n-photon Compton scattering process. A Compton scattering process looks like $k^n e^- \to k e^-$.

You can download this file as a jupyter notebook.

using QEDFeynmanDiagrams

We need some of the packages of the QEDjl-project for base functionality and the ScatteringProcess type.

using QEDcore
using QEDprocesses

Let's decide how many photons our electron interacts with:

n = 4;

Now we setup the scattering process accordingly. We consider all spin/polarization combinations of the particles except for the incoming photons, where the polarizations are synced using QEDbase.SyncedPolarization. This emulates all synced photons having the same, but still indefinite, polarization, for example from a laser.

proc = ScatteringProcess(
    (Electron(), ntuple(_ -> Photon(), n)...),     # incoming particles
    (Electron(), Photon()),                        # outgoing particles
    (AllSpin(), ntuple(_ -> SyncedPol(1), n)...),  # incoming particle spin/pols
    (AllSpin(), AllPol()),                         # outgoing particle spin/pols
)

generic QED process
    incoming: electron (all spins), photon (synced polarization 1), photon (synced polarization 1), photon (synced polarization 1), photon (synced polarization 1)
    outgoing: electron (all spins), photon (all polarizations)

The feynman_diagrams function returns an iterator for all possible Feynman diagrams for this scattering process. With its length overload, we can check how many diagrams there are. For an n-photon Compton process with n incoming photons, this should be $(n+1)!$.

length(feynman_diagrams(proc))

120

Next, we can generate the DAG representing the computation for our scattering process' squared matrix element. This uses ComputableDAGs.jl.

dag = graph(proc)

Graph:
  Nodes: Total: 2262, QEDFeynmanDiagrams.ComputeTask_PropagatePairs: 152, QEDFeynmanDiagrams.ComputeTask_CollectTriples: 16, 
         QEDFeynmanDiagrams.ComputeTask_Triple: 480, QEDFeynmanDiagrams.ComputeTask_BaseState: 14, QEDFeynmanDiagrams.ComputeTask_CollectPairs: 152, 
         QEDFeynmanDiagrams.ComputeTask_Propagator: 30, QEDFeynmanDiagrams.ComputeTask_Pair: 264, ComputableDAGs.DataTask: 1153, 
         QEDFeynmanDiagrams.ComputeTask_SpinPolCumulation: 1
  Edges: 4185
  Total Compute Effort: 0.0
  Total Data Transfer: 0.0
  Total Compute Intensity: 0.0

In this graph output you can see the number of nodes necessary to compute. Note that for larger processes, the number of total nodes can be lower than the number of Feynman diagrams, even with the added complexity of considering multiple spin and polarization combinations. This is the result of efficient reuse of reappearing parts of Feynman diagrams.

To continue, we will need ComputableDAGs.jl. Since ComputableDAGs.jl uses RuntimeGeneratedFunctions as the return type of ComputableDAGs.get_compute_function, we need to initialize it in our current module.

using ComputableDAGs
using RuntimeGeneratedFunctions
RuntimeGeneratedFunctions.init(@__MODULE__)

With the DAG, the process, and RuntimeGeneratedFunctions initalized, we can now generate the actual computable function:

func = get_compute_function(dag, proc, cpu_st(), @__MODULE__);

Now we need an input for the function, which is a QEDcore.PhaseSpacePoint. For now, we generate random momenta for every particle. In the future, QEDevents will be able to generate physical PhaseSpacePoints.

psp = PhaseSpacePoint(
    proc,
    PerturbativeQED(),
    PhasespaceDefinition(SphericalCoordinateSystem(), ElectronRestFrame()),
    tuple((rand(SFourMomentum) for _ in 1:number_incoming_particles(proc))...),
    tuple((rand(SFourMomentum) for _ in 1:number_outgoing_particles(proc))...),
)

PhaseSpacePoint:
    process: generic QED process "ekkkk -> ek"
    model: perturbative QED
    phasespace definition: spherical coordinates in electron rest frame
    incoming particles:
     -> incoming electron: [0.2403111784641755, 0.10282727961854932, 0.3424054542029059, 0.16514997956542432]
     -> incoming photon: [0.24413961136472517, 0.39348607023041604, 0.4939253572778628, 0.06869867423481668]
     -> incoming photon: [0.42162078852178564, 0.907663350334842, 0.27572463565563454, 0.4707156831728114]
     -> incoming photon: [0.46108239663024564, 0.745896286298965, 0.07686409879986877, 0.4673589234636355]
     -> incoming photon: [0.2757170548381528, 0.5413009615124905, 0.03054331436500446, 0.7584731879010913]
    outgoing particles:
     -> outgoing electron: [0.6423587351878165, 0.6943580567950984, 0.1798293649262288, 0.3526784165411546]
     -> outgoing photon: [0.2097689073502046, 0.3006046334932786, 0.06483965279668757, 0.48919337058828327]

Finally, we can test that the function actually runs and computes something by simply calling it on the PhaseSpacePoint:

func(psp)

0.01803108147702798

If we want, we can benchmark the execution speed too:

using BenchmarkTools
@benchmark func($psp)

BenchmarkTools.Trial: 10000 samples with 1 evaluation.
 Range (min … max):  43.211 μs … 95.598 μs  ┊ GC (min … max): 0.00% … 0.00%
 Time  (median):     43.651 μs              ┊ GC (median):    0.00%
 Time  (mean ± σ):   44.194 μs ±  2.878 μs  ┊ GC (mean ± σ):  0.00% ± 0.00%

  ▅█▆▂▂                                                       ▁
  ██████▇▆▆▆▇█▆▆▅▄▄▅▅▄▄▄▃▃▃▄▃▃▄▄▄▃▃▁▄▄▄▅▅▇▇▇▆▆▆▅▆▆▆▆▅▅▄▅▅▄▅▅▅ █
  43.2 μs      Histogram: log(frequency) by time      59.9 μs <

 Memory estimate: 256 bytes, allocs estimate: 2.

This page was generated using Literate.jl.