We present a lattice QCD calculation of the pion scalar form factor and its associated radii, with fully controlled systematic uncertainties. Lattice results are computed on a set of 17 gauge ensembles with $N_f=2+1$ Wilson Clover-improved sea quarks. These ensembles cover four values of the lattice spacing between $a=0.049\mathrm{fm}$ and $a=0.086\mathrm{fm}$, a pion mass range of $130 - 350\mathrm{MeV}$ and various physical volumes. A precise evaluation of the challenging quark-disconnected contributions enables an unprecedented momentum resolution, particularly on large and fine ensembles near the physical quark mass.
To ensure reliable extraction of ground-state matrix elements at both zero and non-zero momentum transfer, we use a broad range of source-sink separations, $1.0\ \mathrm{fm} \lesssim t_{\mathrm{sep}} \lesssim 3.25\ \mathrm{fm}$. For the first time, the scalar radii are obtained using a $z$-expansion parametrization of the $Q^2$-dependence of the form factor rather than relying on a simple linear approximation at low momentum transfer.
The physical extrapolation is guided by the next-to-leading order expressions of the scalar radii from chiral perturbation theory. The pertinent low-energy constants are fitted from lattice data, leading to the first lattice determination of $L_4^r$. Systematic uncertainties related to ground state extraction, form factor parametrization, and the physical extrapolation are rigorously quantified using model averaging based on the Akaike Information Criterion (AIC).