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 QEDcore of the QEDjl-project for base functionality and a process type, for which we can use the Mocks submodule from QEDbase for this tutorial. Downstream, a ScatteringProcess from QEDprocesses.jl could be used, for example.

using QEDcore
using QEDbase.Mocks

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 = QEDbase.Mocks.MockProcessSP(
    (Electron(), Photon()),                        # incoming particles
    (Electron(), ntuple(_ -> Photon(), n)...),     # outgoing particles
    (AllSpin(), AllPol()),                         # incoming particle spin/pols
    (AllSpin(), ntuple(_ -> SyncedPol(1), n)...),  # outgoing particle spin/pols
)

QEDbase.Mocks.MockProcessSP{Tuple{Electron, Photon}, Tuple{Electron, Vararg{Photon, 4}}, Tuple{AllSpin, AllPolarization}, Tuple{AllSpin, Vararg{SyncedPolarization{1}, 4}}}((electron, photon), (electron, photon, photon, photon, photon), (all spins, all polarizations), (all spins, synced polarization 1, synced polarization 1, synced polarization 1, synced polarization 1))

The number_of_diagrams function returns how many diagrams there are for a given process. For an n-photon Compton process with n incoming photons, this should be $(n+1)!$.

number_of_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_CollectTriples: 16, QEDFeynmanDiagrams.ComputeTask_PropagatePairs: 152, 
         QEDFeynmanDiagrams.ComputeTask_CollectPairs: 152, QEDFeynmanDiagrams.ComputeTask_SpinPolCumulation: 1, QEDFeynmanDiagrams.ComputeTask_BaseState: 14, 
         QEDFeynmanDiagrams.ComputeTask_Pair: 264, QEDFeynmanDiagrams.ComputeTask_TripleNegated: 480, DataTask: 1153, 
         QEDFeynmanDiagrams.ComputeTask_Propagator: 30
  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__)

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,
    MockModel(),
    FlatPhaseSpaceLayout(TwoBodyRestSystem()),
    tuple((rand(SFourMomentum) for _ in 1:number_incoming_particles(proc))...),
    tuple((rand(SFourMomentum) for _ in 1:number_outgoing_particles(proc))...),
)

PhaseSpacePoint:
    process: QEDbase.Mocks.MockProcessSP{Tuple{Electron, Photon}, Tuple{Electron, Vararg{Photon, 4}}, Tuple{AllSpin, AllPolarization}, Tuple{AllSpin, Vararg{SyncedPolarization{1}, 4}}}((electron, photon), (electron, photon, photon, photon, photon), (all spins, all polarizations), (all spins, synced polarization 1, synced polarization 1, synced polarization 1, synced polarization 1))
    model: mock model
    phase space layout: FlatPhaseSpaceLayout{TwoBodyTargetSystem{Energy{2}}}(TwoBodyTargetSystem{Energy{2}}(Energy{2}()))
    incoming particles:
     -> incoming electron: [0.33777421325532087, 0.4041838155621277, 0.9073041778928463, 0.3535777796262496]
     -> incoming photon: [0.9058584032319295, 0.807899194534859, 0.6741059848547186, 0.8089565230845425]
    outgoing particles:
     -> outgoing electron: [0.810043217420986, 0.04931756589273928, 0.39921229228756483, 0.24398388293140372]
     -> outgoing photon: [0.9492532689397122, 0.349798548188086, 0.21183294464709468, 0.9836159972790699]
     -> outgoing photon: [0.021143838352322475, 0.9760353545389303, 0.3745520477932829, 0.7122535081366017]
     -> outgoing photon: [0.7504742461910624, 0.707003099721219, 0.7962047256822731, 0.0922807226886101]
     -> outgoing photon: [0.9348896798078674, 0.05037240848335256, 0.3751637481379908, 0.28684767681448653]

With the DAG, the process, RuntimeGeneratedFunctions initialized, and an input type to use, we can now generate the actual computable function:

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

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

func(psp)

8.413186146440346

We can benchmark the execution speed too:

using BenchmarkTools
@benchmark func($psp)

BenchmarkTools.Trial: 10000 samples with 1 evaluation per sample.
 Range (min … max):  34.445 μs … 92.323 μs  ┊ GC (min … max): 0.00% … 0.00%
 Time  (median):     34.866 μs              ┊ GC (median):    0.00%
 Time  (mean ± σ):   35.456 μs ±  2.658 μs  ┊ GC (mean ± σ):  0.00% ± 0.00%

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

 Memory estimate: 256 bytes, allocs estimate: 2.

This page was generated using Literate.jl.