  Results are presented for the nucleon isovector electromagnetic form factors
using 11 ensembles generated by the MILC collaboration using the 2+1+1-flavors
HISQ action. They span 4 lattice spacings $a \sim$ 0.06, 0.09, 0.12 and 0.15~fm
and 3 values of $M_\pi \sim 135, 225$ and 315 MeV. High-statistics estimates
are used to perform a simultaneous extrapolation in the lattice spacing,
lattice volume and light-quark masses. The $Q^2$ dependence over the range
0.05-1.4 ${\rm GeV}^2$ is investigated using both the $z$-expansion and the
dipole form. Final $z$-expansion estimates for the isovector r.m.s. radius are
$r_E = 0.769(27)(30)$ fm $r_M = 0.671(48)(76)$ fm and $\mu^{p-n} =
3.939(86)(138)$ Bohr magneton. The first error is the combined uncertainty from
the leading-order analysis, and the second is an estimate of the additional
uncertainty due to using the leading order chiral-continuum-finite-volume fits.
The dipole estimates, $r_E = 0.765(11)(8)$ fm, $r_M = 0.704(21)(29)$ fm and
$\mu^{p-n} = 3.975(84)(125)$, are consistent with those from the $z$-expansion
but with smaller errors. Our analysis highlights three points. First, all data
from the eleven ensembles and existing lattice data on, or close to, physical
mass ensembles from other collaborations collapses more clearly onto a single
curve when plotted versus $Q^2/M_N^2$ as compared to $Q^2$ with the scale set
by quantities other than $M_N$. The difference between these two analyses is
indicative of discretization errors, some of which presumably cancel when the
data are plotted versus $Q^2/M_N^2$. Second, the size of the remaining
deviation of this common curve from the Kelly curve is small and can be
accounted for by statistical and possible systematic uncertainties. Third, to
improve lattice estimates, high statistics data for $Q^2 < 0.1$ ${\rm GeV}^2$
are needed.
