[an error occurred while processing this directive]
Package-X is hosted by Hepforge, IPPP Durham

The Scattering of Light By Light

This tutorial will illustrate a few strategies you can employ to streamline complex calculations with Package-X by using the scattering of light by light in quantum electrodynamics as an example. We will aim to calculate the scattering amplitude and differential cross section. In certain kinematic limits, compact analytic expressions exist and can be quickly obtained. Otherwise, the general result with full kinematic dependence is extremely complicated, and can take a long time to calculate. But with clever strategies it can be obtained in a reasonable amount of time, and can even be numerically evaluated. We will also save the final result so it can be reloaded upon restarting Mathematica.

Setup

At leading order, the Feynman diagrams are:

    

"s-channel""t-channel""u-channel"

In all diagrams, the internal fermions are electron lines, with their momenta directed with the flow of charge. The QED vertex factor is:

By charge conjugation invariance of QED, each diagram in the second row contributes the same as the diagram above it. Therefore, we only need to add the first three diagrams and multiply the result by 2. We will parametrize the amplitude as:

where the individual contributions are

subject to the on shell kinematic conditions , , and .
The external photon polarization vectors which contract into the loop integrals satisfy transversality . Therefore, it is safe to set in the result.

Define on shell conditions and photon polarization transversality conditions:
In[1]:=
Click for copyable input
Out[1]=
In[2]:=
Click for copyable input
Out[2]=
Initiate the integrals (omitting the factor ), and apply the on shell conditions:
In[3]:=
Click for copyable input

Properties of the photon-photon scattering amplitude

The one loop photon-photon scattering amplitude satisfies a number of properties which can be used to check the calculation.

Ultraviolet behavior

The four photon Green function is expected to be ultraviolet finite.

Apply LoopRefine with the option PartUVDivergent to compute only the ultraviolet divergent part of the total amplitude. Considered separately, each diagram is divergent.
In[7]:=
Click for copyable input
Out[7]=
Out[8]=
Out[9]=
The total amplitude is UV finite:
In[10]:=
Click for copyable input
Out[10]=

QED Ward identities

The scattering amplitude is expected to satisfy the on shell Ward identity
at each of the four vertices

Contract the result with tensors p1μ, p2ν, p3ρ, and p4σ, apply the on shell conditions and photon transversality condition apply LoopRefine (about a minute runtime to verify all four identities on a modern machine):
In[11]:=
Click for copyable input
Out[11]=
Out[12]=
Out[13]=
Out[14]=
It is more efficient to verify these by directly contracting the integrands before calling LoopIntegrate (checking p1μ Ward identity):
In[15]:=
Click for copyable input
Out[18]=

Ward identities in Mathematica 10 and above

Mathematica 10 features a new evaluation control function Inactive which stops functions from evaluating until they are enabled with Activate. We can use this feature to check the Ward identities with ease:

In[19]:=
Click for copyable input
In[22]:=
Click for copyable input
In[23]:=
Click for copyable input
In[24]:=
Click for copyable input
In[25]:=
Click for copyable input

Scattering amplitude in various limits

With Package-X, we can compute the light by light scattering amplitude in different kinematic limits. Remember to restore the omitted factor , and the factor implicit in the normalization of LoopIntegrate.

Low energy (Euler-Heisenberg) limit

In the low energy limit , the amplitude starts at :

Use LoopRefineSeries to obtain the large electron mass behavior of the amplitude (about 15 seconds runtime on a modern machine):
In[19]:=
Click for copyable input
Out[19]=

High energy/large momentum transfer limit

In the high energy limit , the electron mass is negligible.

Make the replacement m0 and call LoopRefine (about 20 seconds runtime on a modern machine):
In[20]:=
Click for copyable input
Out[20]=
In the massless limit, individual Passarino-Veltman funcions contributing to the amplitude are infrared divergent. However the divergences cancel in the total, can can be checked by the absence of 1/ϵ poles:
In[21]:=
Click for copyable input
Out[21]=

Forward scattering limit

In the forward limit and so that . By photon polarization transversality, and can be set to zero.

Make the replacement p3p1, p4p2, t0, p1ρ0 and p2σ0. Call LoopRefine to calculate (about 15 seconds runtime on a modern machine):
In[22]:=
Click for copyable input
Out[22]=

General case

Apply LoopRefine on the full amplitude (about 40 seconds runtime on a modern machine):
In[23]:=
Click for copyable input
Out[23]=

Unpolarized differential cross section

The tensor algebraic capabilities of Package-X can also aid in the evaluation of the unpolarized differential cross section. Below, the cross section will be calculated in three kinematic cases: the low energy limit, the high energy limit, and the general case.

The unpolarized photon-photon differential cross section is given by:

where and .

To obtain the matrix element squared , the complex conjugate of the matrix element is needed,

After multiplying, average over the initial and sum over final photon polarizations using

valid for amplitudes satisfying the Ward identity. Then, the matrix element squared is:

Low energy (Euler-Heisenberg) limit

Since the scattering amplitude in the low energy limit is purely real (below all physical threshold), no complex conjugation is needed. The spin-averaged matrix element squared and therefore the differential cross section are readily obtained.

Contract the product of the low energy amplitude with itself, and apply the on shell kinematic conditions. The factor has been restored:
In[24]:=
Click for copyable input
Out[24]=
Out[25]=
Multiply by the flux and phase space factors to obtain the differential cross section:
In[26]:=
Click for copyable input
Out[26]=
In[27]:=
Click for copyable input
Out[27]=

High energy/large momentum transfer limit

In the high energy limit, the amplitude is complex, with the imaginary parts arising from the various logarithmic factors. After selectively extracting the imaginary parts, the conjugate amplitude is constructed from which the cross section follows.

Use Cases and DeleteDuplicates to construct a list of the logarithmic factors appearing in the high energy scattering amplitude:
In[28]:=
Click for copyable input
Out[28]=
In the physical region (, and ), only the logs give an imaginary part. Define the amplitude and its conjugate by explictly separating the imaginary part from the logarithms:
In[29]:=
Click for copyable input
Contract the product of the low energy amplitude with itself, and apply the on shell kinematic conditions. Collect the result by the various logarithmic factors, and Simplify each coefficient.
In[31]:=
Click for copyable input
Out[31]=
Multiply by the flux and phase space factors to obtain the differential cross section:
In[32]:=
Click for copyable input
Out[32]=
The result is not particularly compact, but it can be plotted as a function of the center-of-mass scattering angle:
In[33]:=
Click for copyable input

General case

In the general case, the amplitude is complex, with the imaginary parts arising from the special functions DiscB, ScalarC0, and ScalarD0. The complex conjugate of the scattering amplitude is constructed by selectively conjugating these special functions.

Construct the conjugate amplitude by selectively complex conjugating the special functions in fullAmplitude:
In[34]:=
Click for copyable input
Contract the product of the low energy amplitude with itself, and apply the on shell kinematic conditions (about 10 seconds runtime on a modern machine).
In[35]:=
Click for copyable input
Multiply by the flux and phase space factors to obtain the differential cross section:
In[36]:=
Click for copyable input
The result is extremely complicated, but it can be numerically plotted as a function of the center-of-mass scattering angle (around 80 seconds runtime):
In[37]:=
Click for copyable input
Out[37]=

The reason for the long evaluation time is due to numerous calls to the complicated special functions ScalarC0, and especially ScalarD0 for each kinematic point. The evaluation time can be improved if the function is organized so that each special function appears once in the expression.

There are 5562 instances of ScalarC0 and ScalarD0 each in ampSqFull:
In[38]:=
Click for copyable input
Out[38]=
Out[39]=
Collect the squared amplitude by the special functions, and simplify each coefficient (25 seconds runtime). There are now only 33 instances of the special functions:
In[40]:=
Click for copyable input
Out[41]=
Out[42]=
Construct the differential cross section, and plot it (1 second runtime):
In[43]:=
Click for copyable input
In[44]:=
Click for copyable input
Out[44]=

To avoid having to recalculate and simplify the differential cross section each time Mathematica is restarted, it is useful to save it to disk with DumpSave. The definition can then be loaded with Get when needed across separate sessions.

Save the definition of dσdΩFull in a binary file:
In[52]:=
Click for copyable input
Load the definitions in "QED-LBL-ampSq.mx":
In[54]:=
Click for copyable input