  We present a lattice QCD calculation of the double-virtual neutral pion
transition form factor, with the goal to cover the kinematic range relevant to
hadronic light-by-light scattering in the muon $g-2$. Several improvements have
been made compared to our previous work. First, we take into account the
effects of the strange quark by using the $N_f=2+1$ CLS gauge ensembles.
Secondly, we have implemented the on-shell $\mathcal{O}(a)$-improvement of the
vector current to reduce the discretization effects associated with Wilson
quarks. Finally, in order to have access to a wider range of photon
virtualities, we have computed the transition form factor in a moving frame as
well as in the pion rest-frame. After extrapolating the form factor to the
continuum and to physical quark masses, we compare our results with
phenomenology. We extract the normalization of the form factor with a precision
of 3.5\% and confirm within our uncertainty previous somewhat conflicting
estimates for a low-energy constant that appears in chiral perturbation theory
for the decay $\pi^0 \to \gamma\gamma$ at NLO. With additional input from
experiment and theory, we reproduce recent estimates for the decay width
$\Gamma(\pi^0 \to \gamma\gamma)$. We also study the asymptotic large-$Q^2$
behavior of the transition form factor in the double-virtual case. Finally, we
provide as our main result a more precise model-independent lattice estimate of
the pion-pole contribution to hadronic light-by-light scattering in the muon
$g-2$: $a_{\mu}^{\mathrm{HLbL}; \pi^0} = (59.7 \pm 3.6) \times 10^{-11}$. Using
in addition the normalization of the form factor obtained by the PrimEx
experiment, we get the lattice and data-driven estimate
$a_{\mu}^{\mathrm{HLbL}; \pi^0} = (62.3 \pm 2.3) \times 10^{-11}$.
