  We complement studies of the neutral pion transition form factor pi^0 -->
gamma^(*) gamma^(*) with calculations for the electromagnetic decay widths of
the processes pi^0 --> e^+ e^-, pi^0 --> e^+ e^- gamma and pi^0 --> e^+ e^- e^+
e^-. Their common feature is that the singly- or doubly-virtual transition form
factor serves as a vital input that is tested in the non-perturbative
low-momentum region of QCD. We determine this form factor from a
well-established and symmetry-preserving truncation of the Dyson-Schwinger
equations. Our results for the three- and four-body decays match results of
previous theoretical calculations and experimental measurements. For the rare
decay we employ a numerical method to calculate the process directly by
deforming integration contours, which in principle can be generalized to
arbitrary integrals as long as the analytic structure of the integrands are
known. Our result for the rare decay is in agreement with dispersive
calculations but still leaves a 2 sigma discrepancy between theory and
experiment.
