feynsage: a Tutorial for Someone Starting QFT

Drawing made by me.

In a first QFT course the calculations arrive in a fixed order. First you learn to take traces of gamma matrices. Then you square an amplitude, average over spins, turn it into a cross section you could compare with an experiment. Loops come later, once the easy things are safe.

feynsage is a SageMath package that does these calculations. This page follows that same order.

Who this is for. You have met the Dirac equation, you have seen a Feynman diagram, you know that a cross section is what an experiment measures. You do not need to have taken a single trace by hand, because we start there.

Everything on this page is real output. Each block was run, the answer captured, then typeset.

Installing it

git clone https://github.com/aburousan/feynsage.git
cd feynsage
./install.sh          # add --test to run the test suite afterwards

The script finds SageMath with FORM, or installs them. Then, inside Sage:

from feynsage import *

What else is in there

This page uses a small corner of the package. The functions fall into a few groups, so you know what to look for when you need it.

info(f) prints what any one of them computes with what its output means, so you can explore without the reference open.

Part 1: gamma matrices and traces

What these objects are

The Dirac equation describes a spin one-half particle such as an electron. To write it down you need four matrices, \(\gamma^0,\gamma^1,\gamma^2,\gamma^3\), each one a \(4\times4\) matrix. They are called gamma matrices or Dirac matrices.

You almost never need the matrices themselves. Everything follows from one rule, the anticommutation relation

Two more pieces of notation. A slash is a momentum contracted with the gammas,

and the trace of a matrix is the sum of its diagonal entries, written \(\mathrm{tr}\).

Taking traces

In feynsage you declare your momenta with your indices, then write the trace the way the book writes it. FORM does the algebra underneath, you never see it.

p, q = momenta('p q')
mu, nu, rho, sigma_ = lorentz_indices('mu nu rho sigma_')

Start with the simplest one, the trace of the identity matrix. It is a \(4\times4\) matrix, so:

dirac_trace(one())

Two gammas:

dirac_trace(gamma(mu)*gamma(nu))

An odd number of them gives nothing at all:

dirac_trace(gamma(mu)*gamma(nu)*gamma(rho))

Four gammas, the identity that takes half a page by hand:

dirac_trace(gamma(mu)*gamma(nu)*gamma(rho)*gamma(sigma_))

Those four results are exactly the trace theorems Peskin and Schroeder collect in equation (5.5) on page 134, which is the box a first course asks you to memorise. Now with momenta in them:

dirac_trace(slash(p)*slash(q))

dirac_trace(slash(p)*gamma(mu)*slash(q)*gamma(nu))

Each of these came back in under three hundredths of a second. The last one is the trace you meet in the first real calculation of the course, which is where we go next.

Part 2: a cross section you could measure

The process

Take the reaction an electron-positron collider actually does:

An electron meets a positron, they annihilate into a photon, the photon turns into a muon with an antimuon. It is the standard first calculation because it has exactly one diagram in QED, no complications.

Set the process up. We ask for pure QED, so only the photon is exchanged, with the electron and muon masses set to zero since the collider energy is far above both:

qed = QED()
P = process('e- e+ -> mu- mu+', model=qed, massless=['e', 'mu'])
len(P.diagrams)

One diagram, as expected.

The squared amplitude

P.squared()

That is \(|\mathcal M|^2\), already averaged over the two initial spins, summed over the final ones. Written tidily,

This is Peskin's equation (5.10) on page 135, taken with the muon mass set to zero. feynsage built the diagram, put in the spin sums, took the traces of Part 1, contracted the indices, then averaged, in two hundredths of a second.

Where the shape comes from: one spin state at a time

Before the cross section, one more thing feynsage can do that makes the answer less mysterious.

The \(|\mathcal M|^2\) above was averaged over spins, so the spin information is gone. Ask instead for a single spin state. For a massless fermion the spin along the direction of motion is called helicity, either \(+1\) or \(-1\):

for he in (1, -1):
    for hm in (1, -1):
        P.squared(helicities={'e-': he, 'e+': -he, 'mu-': hm, 'mu+': -hm}).factor()

These are Peskin's equations (5.21) to (5.24) on page 143.

Look at what the table says. Each single helicity state gives a lopsided answer, \((1+\cos\theta)^2\) throwing everything forwards or \((1-\cos\theta)^2\) throwing it backwards. Neither is symmetric. Only when you average the four of them,

does the familiar shape appear. The \(1+\cos^2\theta\) you plot is the sum of two lopsided pieces, not a basic fact about the process.

P.squared(helicities={'e-': +1, 'e+': +1, 'mu-': +1, 'mu+': -1}).factor()

Turning it into a cross section

A squared amplitude is not yet something you can measure. You must divide by the flux, integrate over the directions the muons can come out in. The useful quantity is the differential cross section \(\mathrm d\sigma/\mathrm d\cos\theta\), which says how much of the scattering goes into each angle, with \(\theta\) the angle between the incoming electron with the outgoing muon.

P.dsigma_dcos()

Read the shape off it: the answer is proportional to \(1+\cos^2\theta\). Using \(e^2 = 4\pi\alpha\) this is

This is Peskin's equation (5.14) on page 137, written there as \(\mathrm d\sigma/\mathrm d\Omega = (\alpha^2/4E_{\rm cm}^2)(1+\cos^2\theta)\). The two agree: the solid angle brings a factor \(2\pi\), since nothing depends on the azimuth, with \(E_{\rm cm}^2 = s\).

More muons come out along the beam than across it, in that particular shape. It is a real prediction that experiments have confirmed many times.

The angular distribution at \(\sqrt{s}=10\) GeV, in picobarns. The solid curve is what feynsage computed, the dashed one is the shape \(1+\cos^2\theta\) normalised at the middle. They lie on top of each other.

Now integrate over all angles to get one number:

P.sigma(sqrt_s=10, unit='pb')

Check it against the textbook. In the high-energy limit the total cross section is

which is the formula on Peskin page 137, where the text notes that dimensional analysis fixes everything except the factor \(4\pi/3\). Put the numbers in by hand:

float(4*pi*(1/137.035999084)^2/(3*10^2)*0.3893793721e9)

against feynsage's \(868.5447687569202\). The two agree in all sixteen digits but the last, which is rounding. The \(0.3893793721\times10^{9}\) is the number of picobarns in one \(\mathrm{GeV}^{-2}\).

That is a complete calculation, from naming a reaction to a number with units, in four lines.

The diagram itself

feynsage draws what it built, so you can check you set up the process you meant:

from feynsage.diagram_plot import draw_diagrams
draw_diagrams(P.diagrams, ncols=1, size=2.6)

The electron with the positron meet, make a photon, the photon makes the muon pair. Arrows run along the fermion lines, the wavy line is the photon.

Part 3: the same process at one loop

Everything so far was tree level, which means no closed loops. The next order in perturbation theory adds one loop to every diagram. Ask feynsage for them:

L = P.loop_diagrams()
len(L)

draw_diagrams(L, ncols=4, size=2.3)

All seven one-loop diagrams of this process in QED, drawn by feynsage. (1) and (2) are the boxes, one crossed. (3) corrects the electron vertex, (4) the muon vertex. (5), (6) and (7) put a fermion loop inside the photon line, with an electron, a muon then a tau going round.

The picture already tells you the structure of the calculation. The seven split into three kinds:

• Vertex corrections, (3) and (4): the photon is emitted after the incoming pair has already exchanged one, which changes the strength of the coupling at the vertex.

• Vacuum polarisation, (5) to (7): the photon spends part of its life as a charged pair. This is the one that makes the electric charge depend on energy.

• Boxes, (1) and (2): a photon crosses between the electron side with the muon side, so the two halves are no longer independent.

The one-loop correction

The quantity that corrects the cross section is the interference between the tree amplitude with the one-loop one. virtual() computes it:

V = P.virtual()

It comes back as a Laurent series in \(\epsilon\). Take it relative to the tree result, at ninety degrees, where \(t = u = -s/2\):

ratio = (V/P.squared()).subs({t: -s/2, u: -s/2}).simplify_full()
pole_parts(ratio)

The coefficient of \(1/\epsilon^2\) is

which is \(-\alpha/\pi\), since \(e^2 = 4\pi\alpha\). The coefficient of \(1/\epsilon\) is

What survives: the charge depends on energy

Of the three kinds of diagram, the vacuum polarisation leaves a mark you can measure. Isolate it by computing the photon self-energy on its own:

qed = QED(muons=False, taus=False)
p = momenta('pp_'); P2 = var('P2_')
se = topologies(1, 1, 1, exclude=('tadpoles', 'wf'))
dA = insert_fields(se, ['A'], ['A'], qed)
gM = sum(dA.loop_amplitude([p], [p],
         lambda n, e: n*T.metric(e[0], e[1]), {'pp_^2': 'P2_'}))
Pi2 = gM/((3 - 2*eps)*P2)
uv_part(Pi2).simplify_full()

So the photon self-energy carries the ultraviolet pole

This is Peskin's equation (7.90) on page 252, converted to the \(d = 4-2\epsilon\) convention. Peskin defines \(\epsilon = 4-d\), so his pole is twice this one. The physics is identical, the letter is not.

That single number is why the fine structure constant is not constant. Absorbing the pole into the charge leaves behind a logarithm of the energy, so the charge you measure grows as you probe harder: \(1/137\) at low energy, about \(1/128\) at the mass of the \(Z\). The whole running of the coupling in QED comes out of this one coefficient.

Part 4: loop integrals by themselves

Behind every diagram in Part 3 sits an integral over the momentum running round the loop, which is not fixed by the incoming particles:

These integrals are the hard part of QFT. Many of them come out infinite, which is what renormalisation exists to handle.

feynsage knows the standard ones by name. The simplest, one propagator in a circle, is called the tadpole:

A0('m')

The \(1/\epsilon\) is the infinity, held in a form you can cancel later. \(\mu\) is a reference energy that the method needs. With two propagators you get the bubble, which depends on the momentum \(s\) flowing through:

s = var('s')
B0(s, 'm', 'm')

\(\Lambda\), printed as DiscB, holds all the dependence on \(s\).

Why the bubble has an imaginary part

Here is the one loop result worth seeing even on a first pass. Take the finite part of the bubble, evaluate it along the real \(s\) axis:

The finite part of \(B_0(s;m,m)\) with \(m=\mu=1\). Below \(s=4m^2\) it is real. At \(s=4m^2\) an imaginary part switches on, which is exactly the energy needed to make the two particles inside the loop real instead of virtual. The dashed line is \(\pi\sqrt{1-4m^2/s}\), sitting on what feynsage computed.

The threshold is the optical theorem: the imaginary part of a loop is the rate for actually producing what is inside it, so it must be zero until there is enough energy. Checked at four points:

The conventions, before you compare with a book

Where to go next

feynsage goes a long way past this page: tensor reduction of loop integrals, integral families, IBP reduction with symmetries found by itself, master integrals, Mellin-Barnes, sector decomposition, differential equations, diagrams generated from a model the way FeynArts does it.

• The repository with the full reference: github.com/aburousan/feynsage

• A longer tutorial ships inside it at docs/tutorial/feynsage_tutorial.pdf

• docs/dirac_algebra.md if traces are what you came for

• Many functions explain themselves. Pass explain=True to the one-loop functions such as A0, or call info(f) on a function, to have it print what the output means.

One honest word on scale. feynsage is built for small and medium problems, where every step stays exact, can be checked, explains itself. For very large reductions with millions of equations, Kira with FIRE are the right tools.