pion charge radius defined by $$\langle r^2\rangle =6\frac{dG_\pi(q^2)}{dq^2}\vert_{q^2=0}.$$ It has been known for some time from chiral perturbation theory analysis [@gasserleutwyler; @gasserleutwyler_2] that this quantity diverges logarithmically with vanishing pion mass squared. Quantitative confirmation of such a behavior would provide an important check on the control of chiral behavior in lattice QCD simulations toward the physical point.

The pioneering lattice QCD calculations of the pion form factor appeared more than 20 years ago [@martinelli; @draper], and a number of studies were carried out over the years. Recently, with the development of simulations with dynamical quarks, several groups have attempted calculations with $N_f=2$ [@brommel; @ETMC; @JLQCD] and $N_f=2+1$ [@RBC] dynamical flavors using a variety of quark actions. The $N_f=2$ calculations employed Wilson-clover [@brommel], twisted mass [@ETMC] or overlap [@JLQCD] quark action, and explored the pion mass region down to $M_\pi\approx 300$ MeV. The values for $\langle r^2\rangle$ from those simulations are significantly smaller than the experimental value, and NNLO fits of SU(2) chiral perturbation theory were needed to find consistency with it at the physical pion mass. For $N_f=2+1$ dynamical flavors, there has been a single calculation employing domain-wall quark action [@RBC], which made measurements at a single pion mass of 