We present results for the isovector (p−n) electromagnetic form factors of the nucleon using eleven ensembles of gauge configurations generated by the MILC collaboration using the highly improved staggered quark (HISQ) action with 2+1+1 dynamical flavors. These ensembles span four lattice spacings a ≈ 0.06, 0.09, 0.12 and 0.15 fm and three values of the light-quark masses corresponding to the pion masses Mπ ≈ 135, 225 and 315 MeV. High-statistics estimates using the truncated solver method method allow us to quantify various systematic uncertainties and perform a simultaneous extrapolation in the lattice spacing, lattice volume and light-quark masses. We analyze the Q 2 dependence of the form factors calculated over the range 0.05 Q 2 ∼ 1.4 GeV 2 using both the model independent z-expansion and the dipole ansatz. Our final estimates, using the z-expansion fit, for the isovector root-mean-square radius of nucleon are rE = 0.769 (27)(30) fm, rM = 0.671(48)(76) fm and µ 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 estimates from the dipole ansatz, rE = 0.765(11)(8) fm, rM = 0.704(21)(29) fm and µ p−n = 3.975(84)(125) Bohr magneton, are consistent with those from the z-expansion but with smaller errors. Our analysis highlights three points. First, all our data for form factors 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 2 N as compared to Q 2 with the scale set by quantities other than MN . The difference between these two ways of analyzing the data is indicative of discretization errors, some of which presumably cancel when the data are plotted versus Q 2 /M 2 N . 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 for r 2 E , r 2 M and µ, high statistics data for Q 2 < 0.1 GeV 2 are needed.