The current 3$sigma$ discrepancy between experiment and Standard Model predictions for $pi^0 to e^+e^-$ is reconsidered using the Pade Theory for bivariate functions, the Canterbury approximants. This method provides a model-independent data-driven a
pproximation 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 $2sigma$, 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.
