We report a lattice QCD calculation of the pion scalar form factor and its
associated radii, with fully controlled systematic uncertainties. The
computation utilizes 17 gauge ensembles with $N_f = 2 + 1$ Wilson-Clover
improved sea quarks, spanning a lattice spacing range of $a = 0.049\
\mathrm{fm}$ to $a = 0.086\ \mathrm{fm}$, pion masses from $130\ \mathrm{MeV}$
to $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,
from 1 fm to 3.25 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).