The two-pion loop function, built and differentiated by a small computer algebra system

Adam Zettel · September 25, 2026 · written with an agent

L (q) = \frac{w}{q} ln \frac{w + q}{2 m _{π}}, w = 4 m_{π}^{2} + q^{2}; L (q) = 1 + \frac{1}{2} \int_{0}^{1} ln (1 + x (1 - x) \frac{q ^{2}}{m _{π}^{2}}) d x

The loop function $L (q)$ sets the momentum dependence of the two-pion-exchange nuclear force at next-to-leading order. The program writes it as an expression tape, evaluates it, and differentiates the tape symbolically into a second tape. Both are checked against a 30-digit Feynman-parameter integral that never uses the closed form: $L (200 MeV) = 1.14642288565$ and $d L / d q = 1.24075106 \times 1 0^{- 3}$ /MeV there, to $1 0^{- 14}$ and $2 \times 1 0^{- 17}$ /MeV.

Open in the workbench Download .vivOutputs below are from the site build’s run; Run all runs it again in your browser.

Look inside this program

Atlas
The program’s map: every function and which calls which.
There’s no recorded run of this program yet.

In chiral effective field theory the force between two nucleons is a sum of pion exchanges and short-range contact terms. At next-to-leading order the exchange of two pions enters through one function of the momentum transfer $q$ ,

L (q) = \frac{w}{q} ln \frac{w + q}{2 m _{π}}, w = 4 m_{π}^{2} + q^{2},

which multiplies the central, tensor and spin-spin parts of the potential (Machleidt and Entem, Phys. Rep. 503, 1 (2011), Eqs. (4.9)-(4.12), with their average pion mass $m_{π} = 138.039$ MeV). $L (0) = 1$ , and $L$ grows like $ln q$ at large $q$ . As a function of $q^{2}$ it has a branch cut that starts at $q^{2} = - 4 m_{π}^{2}$ , the two-pion threshold, so in coordinate space the force it carries has range $1/ (2 m_{π})$ , half the one-pion range.

Higher orders need many such functions and their derivatives, which is the job of a computer algebra system. This post builds a small one: an expression tape, a list of operations each of which refers to earlier entries, with an evaluator and a differentiator that turns a tape into the tape of its derivative.

The reference, posts/refs/the-pion-loop-function-by-computer-algebra.py, does not use the closed form. It computes the loop from its Feynman-parameter integral (the equation under the title) with mpmath's quadrature at 30 digits, and the derivative by mpmath's numerical differentiation of that integral.

The tape of $L (q)$

The program writes $L$ the way it is printed: $w$ as the square root of $4 m_{π}^{2} + q^{2}$ , then $(w / q) ln ((w + q) / (2 m_{π}))$ . That is 11 entries. Differentiating appends the derivative of every entry, including the derivatives of constants, which are zero; the tape grows to 36 entries. Nothing simplifies the result, so the derivative tape does more arithmetic than a hand-written derivative would, but it is exact.

Printed

n_L: 11

n_all: 36

Tests 1 of 1 pass

the tape of L has 11 entries, and 36 once differentiated (passes)

Values against the Feynman-parameter integral

The reference evaluates $L$ at eight momenta from its integral representation, which comes from the loop diagram itself, not from the closed form. The tape agrees at every one to $1 0^{- 14}$ , from $q = 1$ MeV, where $L - 1 = 4.4 \times 1 0^{- 6}$ , to 1 GeV. At small $q$ the logarithm in the integral expands, and with $t = q^{2} / m_{π}^{2}$ and $\int_{0}^{1} x (1 - x) d x = 1/6$ , $\int_{0}^{1} x^{2} (1 - x)^{2} d x = 1/30$ ,

L = 1 + \frac{t}{12} - \frac{t ^{2}}{120} + O (t^{3}) .

At $q = 1$ MeV, $t^{2} /120 = 2.3 \times 1 0^{- 11}$ , and the tape matches the two terms to $1 0^{- 13}$ .

A third representation is the dispersion integral: $L$ has an imaginary part above the two-pion threshold $μ = 2 m_{π}$ , and

L (q) = 1 + q^{2} \int_{2 m_{π}}^{\infty} \frac{μ ^{2} - 4 m _{π}^{2}}{μ ^{2} ( μ ^{2} + q ^{2} )} d μ .

The program maps $[2 m_{π}, \infty)$ to $[0, 1)$ with $μ = 2 m_{π} / (1 - u)$ and integrates by Simpson's rule on 2000 intervals. The integrand has a square-root zero at threshold, where Simpson's rule loses its order, and the result is off by $4.4 \times 1 0^{- 7}$ at $q = 200$ MeV.

Printed

L_tape: [1, 1.01079, 1.0416, 1.14642, 1.27928, 1.54876, 1.78481, 2.07352]

L_d: 1.14642

disp_err: 4.41822e-07

L1_minus_one: 4.37334e-06

t_sq_term: 2.29515e-11

moment_1: 0.166667

moment_2: 0.0333333

Tests 6 of 6 pass

the tape matches the Feynman-parameter integral at all eight momenta to 1e-14 (passes)
at q = 1 MeV it is 1 + t/12 - t^2/120, t = q^2/m_pi^2, to 1e-13 (passes)
the dispersion integral by Simpson's rule gives L(200) to 1e-6 (passes)
the dispersion integral is off by 4.4e-7 at q = 200 MeV (passes)
at q = 1 MeV, L - 1 is 4.4e-6 and t^2/120 is 2.3e-11 (passes)
the integrals of x(1 - x) and x^2 (1 - x)^2 over [0, 1] are 1/6 and 1/30 (passes)

The derivative

The derivative tape, evaluated at three momenta, matches mpmath's numerical derivative of the Feynman-parameter integral to $2 \times 1 0^{- 17}$ per MeV, rounding in a number of about 0.001. A central difference on the tape with a step of 1 MeV is off by up to $5.3 \times 1 0^{- 9}$ , which is the $h^{2} L^{'''} /6$ error of that step; the symbolic derivative has no step to choose.

Printed

dL_tape: [0.000426174, 0.00124075, 0.00127628]

dL_fd: [0.000426169, 0.00124075, 0.00127628]

dL_tape_err: 1.70762e-17

dL_fd_err: 5.29477e-09

L(q), tapeFeynman-parameter integral (mpmath)dL/dq, derivative tape (per GeV) (right axis)mpmath derivative (per GeV) (right axis)

Tests 4 of 4 pass

the symbolic derivative matches mpmath's derivative of the integral to 1e-15 per MeV (passes)
the derivative tape is within 2e-17 per MeV of mpmath at all three momenta (passes)
the central difference is off by up to 5.3e-9 (passes)
a central difference with a 1 MeV step is off by 2e-10 to 6e-9 (passes)

The tape of L(q)

Values against the Feynman-parameter integral

The derivative

The tape of $L (q)$