  The current 3$\sigma$ discrepancy between experiment and Standard Model
predictions for $\pi^0 \to e^+e^-$ is reconsidered using the Pad\'e Theory for
bivariate functions, the Canterbury approximants. This method provides a
model-independent data-driven approximation to the decay as soon as
experimental data for the doubly virtual $\pi^0$ transition form factor are
available. It also implements the correct QCD constraints of the form factor
both at low- and high-energies. We reassess the Standard Model result
including, for the first time, a systematic error. Our result, BR$(\pi^0 \to
e^+e^-)=6.23(5)\times 10^{-8}$, still represents a discrepancy larger than
$2\sigma$, unsurmountable with our present knowledge of the Standard Model, and
would claim New Physics if the experimental result is confirmed by a new
measurement. Our method also provides the adequate tool to extract the doubly
virtual form factor from experimental data in a straightforward manner. This
measurement would further shrink our error and establish once and for all the
New Physics nature of the discrepancy. In addition, we remark the challenge
this discrepancy poses in the evaluation of the hadronic light-by-light
scattering contribution to the $(g-2)_{\mu}$, specially confronted with the
foreseen accuracy of the forthcoming $(g-2)_{\mu}$ experiments.
