text:
# Nucleon Charge and Magnetization Densities from Sachs Form Factors

## I Introduction

A fundamental test of the QCD confinement mechanism, either on the lattice or in models, is the electromagnetic structure of the nucleon. This electromagnetic structure is reflected by the electric and magnetic form factors, $`G_E(Q^2)`$ and $`G_M(Q^2)`$, measured by electron scattering with spacelike invariant momentum transfer $`Q`$. At low $`Q`$ one can interpret these form factors as Fourier transforms of the nucleon charge and magnetization densities , but the relationship between form factor and density is complicated by recoil as $`Q`$ increases. Although models of nucleon structure can often calculate the form factor directly, it is still desirable to relate form factors to spatial densities because our intuition tends to be grounded more firmly in space than momentum transfer. In this paper we fit charge and magnetization densities to recent nucleon form factor data using a prescription that accounts for nucleon recoil and Lorentz contraction and is compatible with perturbative QCD (pQCD) scaling at large $`Q^2`$.

Early experiments with modest $`Q^2`$ suggested that

$$G_{Ep}\frac{G_{Mp}}{\mu _p}\frac{G_{Mn}}{\mu _n}G_D$$
(1)
where $`G_D(Q^2)=(1+Q^2/\mathrm{\Lambda }^2)^2`$ with $`\mathrm{\Lambda }^2=0.71`$ (GeV/$`c`$)2 is known as the dipole form factor . However, the naive Fourier transform of the dipole form factor produces an exponential density with an unphysical cusp at the origin. Similarly, data for $`G_{En}`$ at low $`Q^2`$ can be described by the Galster parametrization

$$G_{En}(Q^2)\mu _nG_D(Q^2)\frac{A\tau }{1+B\tau }$$
(2)
where $`A`$ and $`B`$ are constants and $`\tau =(Q/2m)^2`$, but direct Fourier transform of this form factor also produces a cusp at the origin. Using a relativistic inversion formula that accounts for the Lorentz contraction along the momentum transfer, Licht and Pagnamenta  obtained a reasonable fit to proton form factors using a Gaussian density with a more realistic behavior in the interior. Ji  obtained similar results using a relativistic soliton model and we used this model in Ref.  to fit data for Sachs form factors. These models offer plausible radial densities, but are not compatible with pQCD scaling unless one imposes somewhat awkward restrictions upon the Fourier transform, as discussed in Sec. II C. Fortunately, a variant proposed by Mitra and Kumari  complies with pQCD scaling automatically, without need of such constraints. We use this method, described as relativistic inversion, to extract nucleon charge and magnetization densities from data for Sachs form factors. Our fitting procedure minimizes model dependence by employing linear expansions in complete radial bases, such as Fourier-Bessel or Laguerre-Gaussian expansions, and estimates uncertainties arising from the limitation of experimental data to a finite range of momentum transfer using methods originally developed to analyze electron scattering by nuclei. Such an analysis produces a good fit to form factor data using a radial density whose error band reflects both the statistical quality of the data and its limited coverage of momentum transfer. Differences between densities obtained using several variations of the inversion formula are described as discrete ambiguities and are analyzed in detail herein.

Data for $`G_{Mp}`$ and $`G_{Mn}`$ with $`Q^2>1`$ (GeV/$`c`$)2 show significant departures from the simple dipole parametrization, but the extraction of $`G_{Ep}`$ using the traditional Rosenbluth method  becomes increasingly difficult as $`Q^2`$ increases because the dominance of the $`G_{Mp}`$ contribution to the cross section increases with $`Q^2`$. Consequently, there are large statistical uncertainties in Rosenbluth data for $`G_{Ep}`$ at $`Q^2>1`$ (GeV/$`c`$)2 and the discrepancies between comparable experiments suggests that systematic errors in the Rosenbluth analysis are often underestimated . More recently, recoil polarization has been used to measure the ratio $`g_p=G_{Ep}/G_{Mp}`$ directly, without need of Rosenbluth separation. In this technique, the components of the nucleon polarization $`\stackrel{}{P}^{}`$ after scattering by a polarized electron beam are measured along the momentum transfer direction, denoted by $`\widehat{z}`$, and in the $`\widehat{x}`$ direction transverse to $`\widehat{z}`$ in the scattering plane. The form factor ratio is then obtained using

$$\frac{P_x^{}}{P_z^{}}=\sqrt{\frac{2ϵ}{\tau (1+ϵ)}}g$$
(3)
where $`ϵ=(1+(1+\tau )2\mathrm{tan}^2\theta _e/2)^1`$ is the transverse polarization of the virtual photon for an electron scattering angle $`\theta _e`$. For the proton, both components can be measured simultaneously using a polarimeter in the focal plane of a magnetic spectrometer, thereby minimizing systematic uncertainties due to beam polarization, analyzing power, and kinematic parameters. The systematic uncertainty due to precession of the proton spin in the magnetic spectrometer is usually much smaller than the systematic uncertainties in combining the absolute cross sections obtained with different kinematical conditions and acceptances that are needed for the Rosenbluth method. Recent data using the recoil polarization technique  have shown a dramatic, almost linear, decrease in $`G_{Ep}/G_{Mp}`$ for $`Q^2>1`$ (GeV/$`c`$)2. It was suggested that those results demonstrate that the proton charge is distributed over a larger volume than its magnetization, but radial densities were not obtained. Our analysis confirms that interpretation quantitatively.

Similar techniques can be used to obtain the neutron form factor ratio, $`g_n=G_{En}/G_{Mn}`$, using either target or recoil polarization, but in the absence of a target with free neutrons one must employ quasifree scattering from a neutron bound in a light nucleus. Detection of a recoil neutron with quasifree kinematics and small missing momentum tends to minimize uncertainties due to nuclear structure and final state interactions . Although considerable care is still needed at low $`Q^2`$, polarization methods offer smaller systematic errors and less model dependence than traditional Rosenbluth analyses of elastic scattering or quasifree knockout. We extracted the neutron charge density from recent polarization data for $`g_n`$ for $`Q^2<1.6`$ (GeV/$`c`$)2 using the relativistic inversion method. Although the form factor data remain consistent with the Galster parametrization over this range of momentum transfer, the charge density obtained by relativistic inversion is considerably softer than that from nonrelativistic inversion of the Galster form factor and does not feature a cusp at the origin. Over the next several years, extending of the experimental range of momentum transfer should substantially reduce the uncertainty in the interior density.

The model is presented in Sec. II, the analysis procedures in Sec. III, and principal results in Sec. IV. In Sec. V we compare our results to another analyses, discuss the extrapolation to higher $`Q^2`$ and the role of discrete ambiguities in fitted densities. We also form combinations of neutron and proton charge densities that in the naive quark model represent the distribution of up and down quarks in the proton. Finally, our conclusions are summarized in Sec. VI.

## II Model

### A Sachs Form Factors

Matrix elements of the nucleon electromagnetic current operator $`J^\mu `$ take the form

$$N(p^{},s^{})|J^\mu |N(p,s)=\overline{u}(p^{},s^{})e\mathrm{\Gamma }^\mu u(p,s)$$
(4)
where $`u`$ is a Dirac spinor, $`p,p^{}`$ are initial and final momenta, $`q=pp^{}`$ is the momentum transfer, $`s,s^{}`$ are spin four-vectors, and where the vertex function

$$\mathrm{\Gamma }^\mu =F_1(Q^2)\gamma ^\mu +\kappa F_2(Q^2)\frac{i\sigma ^{\mu \nu }q_\nu }{2m}$$
(5)
features Dirac and Pauli form factors $`F_1`$ and $`F_2`$ that depend upon the nucleon structure. Here $`e`$ is the elementary charge, $`m`$ is the nucleon mass, $`\kappa `$ is the anomalous part of the magnetic moment, and $`\gamma ^\mu `$ and $`\sigma ^{\mu \nu }`$ are the usual Dirac matrices (e.g., ). The interpretation of these form factors appears simplest in the nucleon Breit frame where the energy transfer vanishes. In this frame the nucleon approaches with initial momentum $`\stackrel{}{q}_B/2`$, receives three-momentum transfer $`\stackrel{}{q}_B`$, and leaves with final momentum $`\stackrel{}{q}_B/2`$. Thus, the nucleon Breit frame momentum is defined by $`q_B^2=Q^2=q^2/(1+\tau )`$ where $`(\omega ,\stackrel{}{q})`$ is the momentum transfer in the laboratory, $`Q^2=q^2\omega ^2`$ is the spacelike invariant four-momentum transfer, and $`\tau =Q^2/4m^2`$. In the Breit frame for a particular value of $`Q^2`$, the current separates into electric and magnetic contributions

$$\overline{u}(p^{},s^{})\mathrm{\Gamma }^\mu u(p,s)=\chi _s^{}^{}\left(G_E+\frac{i\stackrel{}{\sigma }\times \stackrel{}{q}_B}{2m}G_M\right)\chi _s$$
(6)
where $`\chi _s`$ is a two-component Pauli spinor and where the Sachs form factors are given by

$`G_E`$ $`=`$ $`F_1\tau \kappa F_2`$ (8)
$`G_M`$ $`=`$ $`F_1+\kappa F_2`$ (9)
The similarity of Eq. (6) to the classical current density

$$J^{NR}=(e\rho _{ch}^{NR},\mu \stackrel{}{\sigma }\times \stackrel{}{}\rho _m^{NR})$$
(10)
suggests an identification of charge and magnetization densities

$`\rho _{ch}^{NR}(r)`$ $`=`$ $`{\displaystyle \frac{2}{\pi }}{\displaystyle _0^{\mathrm{}}}𝑑QQ^2j_0(Qr)G_E(Q^2)`$ (12)
$`\mu \rho _m^{NR}(r)`$ $`=`$ $`{\displaystyle \frac{2}{\pi }}{\displaystyle _0^{\mathrm{}}}𝑑QQ^2j_0(Qr)G_M(Q^2)`$ (13)
where $`\mu =1+\kappa `$ is the appropriate static magnetic moment (either $`\mu _p`$ or $`\mu _n`$) relative to the nuclear magneton. However, this naive inversion procedure is described as nonrelativistic (NR) because it ignores the variation of the Breit frame with $`Q^2`$.

### B Intrinsic Form Factors

Let $`\rho _{ch}(r)`$ and $`\rho _m(r)`$ represent spherical charge and magnetization densities in the nucleon rest frame. It is convenient to normalize these intrinsic densities according to

$`{\displaystyle _0^{\mathrm{}}}𝑑rr^2\rho _{ch}(r)`$ $`=`$ $`Z`$ (15)
$`{\displaystyle _0^{\mathrm{}}}𝑑rr^2\rho _m(r)`$ $`=`$ $`1`$ (16)
where $`Z=0,1`$ is the nucleon charge. Fourier-Bessel transforms of the intrinsic densities are defined by

$$\stackrel{~}{\rho }(k)=_0^{\mathrm{}}𝑑rr^2j_0(kr)\rho (r)$$
(17)
where $`k`$ is the spatial frequency (or wave number). We describe $`\stackrel{~}{\rho }(k)`$ as an intrinsic form factor. If one knew how to obtain $`\stackrel{~}{\rho }(k)`$ from data for the appropriate Sachs form factor, the intrinsic density could be obtained simply by inverting the Fourier transform, such that

$$\rho (r)=\frac{2}{\pi }_0^{\mathrm{}}𝑑kk^2j_0(kr)\stackrel{~}{\rho }(k)$$
(18)
The naive nonrelativistic inversion method assumes that $`kQ`$ and $`\stackrel{~}{\rho }(Q)G(Q^2)`$ where $`G(Q^2)`$ is the appropriate Sachs form factor. However, this inversion procedure produces unsatisfactory results for the common dipole and Galster parametrizations — the corresponding radial densities have unphysical cusps at the origin and rather hard cores. For example, the naive Fourier transform of the dipole form factor produces an exponential density. (Although it appears much more complicated, the Galster density can also be obtained in closed form and displays similarly unrealistic behavior near the origin.) Licht and Pagnamenta  attributed these failures of nonrelativistic inversion to the replacement of the intrinsic spatial frequency $`k`$ with the momentum transfer $`Q`$ and demonstrated that by applying a boost from the Breit frame with momentum $`q_B=Q`$ to the rest frame, inversion of the dipole form factor using a reduced spatial frequency

$$k^2=\frac{Q^2}{1+\tau }$$
(19)
softens the density. In fact, a good fit to the data for $`G_{Ep}`$ could then be obtained using a Gaussian density typical of quark models.

Unfortunately, unique relativistic relationships between the Sachs form factors measured by electron scattering at finite $`Q^2`$ and the static charge and magnetization densities in the nucleon rest frame do not exist. The basic problem is that electron scattering measures transition matrix elements between states of a composite system that have different momenta and the transition densities between such states are different from the static densities in the rest frame. Furthermore, the boost operator for a composite system depends upon the interactions among its constituents. Nevertheless, a wide variety of models have employed similar relativistic prescriptions to relate elastic form factors to ground-state densities. The first proposal was made by Licht and Pagnamenta  using a cluster model and a kinematic boost that neglects interactions. The transition form factors were evaluated using the impulse approximation and neglecting relative motion. Mitra and Kumari  obtained similar results using a kinematic transformation that is more symmetric between initial and final states and can be applied to inelastic scattering also. Ji  also obtained similar results using a relativistic Skyrmion model based upon a Lorentz invariant Lagrangian density for which the classical soliton solution can be evaluated in any frame. Quantum fluctuations were then evaluated after the boost. Although an approximation is still required to evaluate the transition form factors, it was argued that this approximation is best in the Breit frame. Holzwarth  extended the soliton model to the timelike regime and analyzed the superconvergence relations needed to obtain spectral functions.

Each of these prescriptions can be represented in the form

$`\stackrel{~}{\rho }_{ch}(k)`$ $`=`$ $`G_E(Q^2)(1+\tau )^{\lambda _E}`$ (21)
$`\mu \stackrel{~}{\rho }_m(k)`$ $`=`$ $`G_M(Q^2)(1+\tau )^{\lambda _M}`$ (22)
where $`G(Q^2)`$ is one of the four Sachs form factors, $`k`$ is the intrinsic spatial frequency given by Eq. (19), and $`\lambda `$ is a model-dependent constant. The most important relativistic effect is Lorentz contraction of spatial distributions in the Breit frame and the corresponding increase of spatial frequency represented by the factor of $`(1+\tau )`$ in Eq. (19). A measurement with Breit-frame momentum transfer $`q_B=Q`$ probes a reduced spatial frequency $`k`$ in the rest frame. The Sachs form factor for a large invariant momentum transfer $`Q^2`$ is determined by a much smaller spatial frequency $`k^2=Q^2/(1+\tau )`$ and thus declines much less rapidly with respect to $`Q^2`$ than the Fourier transform of the density declines with respect to $`k^2`$. In fact, the accessible spatial frequency is limited to $`k2m`$ such that the asymptotic Sachs form factors in the limit $`Q^2\mathrm{}`$ are determined by the Fourier transform of intrinsic densities in the immediate vicinity of the limiting frequency $`k_m=2m`$. In this model, no information can be obtained beyond the limiting frequency determined by the nucleon Compton wavelength. This limitation can be understood as a consequence of relativistic position fluctuations, known as of zitterbewegung, that smooth out radial variations on scales smaller than the Compton wavelength.

Ji  derived $`\lambda _E=0`$ for electric and $`\lambda _M=1`$ for magnetic form factors in the soliton model and attributed the difference between $`\lambda _E`$ and $`\lambda _M`$ to the Lorentz transformation properties of scalar and vector densities. The same choices were employed by Holzwarth . On the other hand, Licht and Pagnamenta  obtained $`\lambda _E=\lambda _M=1`$ using the cluster model, but Mitra and Kumari  found that a more symmetric treatment of the kinematics gives values $`\lambda _E=\lambda _M=2`$ that automatically satisfy the perturbative QCD scaling relations at very large $`Q^2`$. For most of the present analysis we will employ the symmetric choice $`\lambda _E=\lambda _M=2`$ because fewer restrictions upon the behavior of $`\stackrel{~}{\rho }(k)`$ are needed near $`k_m`$ to ensure compatibility with the asymptotic behavior of Sachs form factors expected from dimensional scaling. We describe the variation of a fitted density with the choice of $`\lambda `$ as a discrete ambiguity. The effect of discrete ambiguities upon fitted densities will be examined in Sec. V D.

### C Asymptotic Behavior

The present form factor model with $`\lambda 0`$ suggests that the asymptotic behavior for $`Q2m`$ is given by

$`G(Q^2)`$ $``$ $`\left({\displaystyle \frac{k_m}{Q}}\right)^{2\lambda }(\stackrel{~}{\rho }(k_m)(2\lambda \stackrel{~}{\rho }(k_m)+k_m\stackrel{~}{\rho }^{}(k_m)){\displaystyle \frac{k_m^2}{2Q^2}}`$ (23)
$`+`$ $`(4\lambda (\lambda +1)\stackrel{~}{\rho }(k_m)+(3+4\lambda )k_m\stackrel{~}{\rho }^{}(k_m)+k_m^2\stackrel{~}{\rho }^{\prime \prime }(k_m)){\displaystyle \frac{k_m^4}{8Q^4}}+\mathrm{})`$ (24)
where $`\stackrel{~}{\rho }^{}`$ and $`\stackrel{~}{\rho }^{\prime \prime }`$ are derivatives of the momentum-space density evaluated at the limiting frequency $`k_m=2m`$. Evidently, noninteger values of $`\lambda `$ are incompatible with the perturbative QCD prediction  that $`GQ^4`$ aside from logarithmic corrections. Similarly, $`\lambda >2`$ is excluded also, leaving just three choices. If we choose $`\lambda =2`$, then we need only require $`\stackrel{~}{\rho }(k_m)0`$ to obtain consistency with pQCD. Thus, the proposal by Mitra and Kumari  of $`\lambda _E=\lambda _M=2`$ offers the most natural approach to the pQCD limit. If we choose $`\lambda =1`$, as recommended by Refs.  for $`G_M`$ or by Ref.  for both $`G_E`$ and $`G_M`$, then we must require $`\stackrel{~}{\rho }(k)`$ to have a node at $`k_m`$ such that $`\stackrel{~}{\rho }(k_m)=0`$ and $`\stackrel{~}{\rho }^{}(k_m)0`$. Finally, if we choose $`\lambda =0`$ as recommended by Refs.  for $`G_E`$, then we must impose the somewhat unnatural constraints $`\stackrel{~}{\rho }(k_m)=\stackrel{~}{\rho }^{}(k_m)=0`$ with $`\stackrel{~}{\rho }^{\prime \prime }(k_m)0`$. Thus, it appears that the usefulness of the chiral soliton model is limited to $`Q^24m^2`$ and in order to fit data for larger $`Q^2`$ with that model Holzwarth found it necessary to artificially increase the soliton mass .

Although the intrinsic form factors $`\stackrel{~}{\rho }(k)`$ obtained using either dipole or Galster functions for $`G(Q^2)`$ in Eq. (II B) are compatible with the pQCD constraints upon $`\stackrel{~}{\rho }(k_m)`$, neither can be inverted using Eq. (18) with $`\lambda =0`$ because $`\stackrel{~}{\rho }(k)G(k_m^2)(k_m/k)^{2\lambda }`$ for $`k\mathrm{}`$. The inversion integrals for these functions converge well for $`\lambda =2`$, slowly for $`\lambda =1`$, and diverge for $`\lambda =0`$. Recognizing that pQCD favors $`\lambda =2`$, we expect $`\stackrel{~}{\rho }(k)`$ to have an asymptotic $`k^4`$ behavior with an amplitude determined by the nucleon-antinucleon annihilation process $`N\overline{N}e^{}e^+`$ at threshold. Similarly, the behavior of $`\stackrel{~}{\rho }(k)`$ for $`k\mathrm{}`$ should be determined by the electromagnetic annihilation data for $`Q^24m^2`$. However, we have not attempted to incorporate electromagnetic data for timelike $`Q^2`$ in the present analysis because it is not clear that the prescription for intrinsic form factors should apply to that regime. A more general analysis of the analytic structure of the form factors can be made using dispersion theory , but that approach does not consider the densities that are the subject of the present analysis. Nor do we consider here the modifications of pQCD scaling due to logarithmic running of the strong coupling .

With the exception of $`G_{Mp}`$, the available data constrain $`\stackrel{~}{\rho }(k)`$ very little near $`k_m`$ because the ratio

$$\frac{k}{k_m}1\frac{k_m^2}{2Q^2}$$
(25)
approaches unity relatively slowly as $`Q^2`$ increases. Thus, we will find that the choice of $`\lambda `$ has very little effect upon the fit to data for Sachs form factors, but does have a strong effect upon the extrapolation of fitted form factors beyond the measured range of $`Q^2`$. By incorporating pQCD scaling in its basic parametrization, the choice $`\lambda =2`$ limits the range of variation available to extrapolated form factors. Conversely, without explicit enforcement of pQCD by means of somewhat artificial constraints upon $`\stackrel{~}{\rho }(k)`$ near $`k_m`$, fits with $`\lambda <2`$ permit much wider latitude at large $`Q^2`$. The data for $`G_{Mp}`$ for $`Q^2>20`$ (GeV/$`c`$)2 exhibit scaling and automatically enforce the appropriate constraints upon $`\stackrel{~}{\rho }(k_m)`$, but data available for the other three electromagnetic form factors do not. Consequently, one could impose constraints upon $`\stackrel{~}{\rho }(k_m)`$ with little effect upon the fits in the measured range of $`Q^2`$. However, we chose not to employ constraints of this kind, which seem rather artificial, and to permit fits with $`\lambda <2`$ the greatest possible latitude.

We use $`\lambda _E=\lambda _M=2`$ for most of the present work, but will discuss the consequences of the discrete ambiguity in Sec. V D. Note that our previous work , motivated by the soliton model, used $`\lambda _E=0`$ and $`\lambda _M=1`$.

### D Moments

It is customary to describe the low $`Q^2`$ behavior of a form factor in terms of a transition radius obtained from integral moments of the underlying density, but care must be taken with the relativistic relationship between a Sachs form factor and its intrinsic density. We define integral moments by

$$M_\alpha =_0^{\mathrm{}}𝑑rr^{2+\alpha }\rho (r)$$
(26)
where $`\alpha `$ is an even integer. For a charge density these moments are related to the electric form factor by

$`M_0`$ $`=`$ $`G(0)`$ (28)
$`M_2`$ $`=`$ $`6{\displaystyle \frac{dG(Q^2)}{dQ^2}}|_{Q^20}{\displaystyle \frac{3\lambda }{2m^2}}G(0)`$ (29)
while for magnetization we divide by the magnetic moment. Thus, one expects $`M_0=Z`$ for charge densities and $`M_0=1`$ for magnetization densities. Notice that the lowest nonvanishing moment is free of discrete ambiguities, but that higher moments depend upon $`\lambda `$. For example, the mean-square neutron charge radius reduces to

$$r^2_n=6\frac{dG(Q^2)}{dQ^2}|_{Q^20}$$
(30)
because the charge vanishes, while the proton radius retains a small dependence upon $`\lambda `$

$$r^2_{\lambda ,p}=r^2_{0,p}\frac{3\lambda }{2m_p^2}$$
(31)
due to the discrete ambiguity in the intrinsic density. This term, equal to $`0.066\lambda `$ in units of fm2, appears to be similar to the famous Foldy contribution to the neutron charge radius  but has a different origin because it does not depend upon the anomalous magnetic moment.

When $`M_00`$ it is useful to distinguish between a radius parameter

$$\xi =\left(6\frac{d\mathrm{ln}G(Q^2)}{dQ^2}\right)_{Q^20}^{1/2}=\left(\frac{M_2}{M_0}+\frac{3\lambda }{2m^2}\right)^{1/2}$$
(32)
based upon the initial logarithmic derivative of a Sachs form factor and the rms radius

$$r^2_\lambda ^{1/2}=\left(\xi ^2\frac{3\lambda }{2m^2}\right)^{1/2}$$
(33)
of the corresponding density obtained for specified $`\lambda `$. Thus, $`\xi `$ is a model-independent property of the form factor data while $`r^2^{1/2}`$ is subject to a discrete ambiguity. These radii agree for $`\lambda =0`$, but $`r^2^{1/2}`$ is smaller than $`\xi `$ for larger $`\lambda `$ due to the zitterbewegung correction.

Accurate calculations for many phenomena in atomic physics, such as the Lamb shift, require corrections for the finite size of nucleons. Although it might appear that the nucleon size should be determined by an integral moment of a nucleon density through Eq. (26), the static radial density is not directly measurable by electron scattering. The discrete ambiguity between the initial slope of the Sachs form factor and its associated transition density reflects the model-dependence of the relativistic inversion procedure arising from the treatment of zitterbewegung. By convention, QED theorists have decided to identify the radius with the initial slope of the Sachs form factor and to treat recoil, vacuum polarization, zitterbewegung, and other effects as separate corrections (e.g. ). To distinguish between various determinations of nucleon size, we describe the model-independent quantity $`\xi `$ as the Sachs radius and the model-dependent rms radius obtained from moments of a fitted radial density as an intrinsic radius. When necessary, these radii are further qualified as either charge or magnetic and for intrinsic radii by the value of $`\lambda `$. The Sachs charge radius is usually the most appropriate for QED applications.

### E Example: Gaussian density

It is instructive to consider the Sachs form factor that would be produced by a simple Gaussian density

$$\rho (r)=\frac{4}{b^3\sqrt{\pi }}\mathrm{exp}((r/b)^2)\stackrel{~}{\rho }(k)=\mathrm{exp}((kb/2)^2)$$
(34)
that is typical of quark models. The form factor obtained using Eq. (II B) with $`\lambda =2`$ is compared with the familiar dipole form factor in Fig. 1 for several choices of $`b`$; note with $`b=0.556`$ fm the Gaussian parametrization has the same rms radius as the dipole form factor. These curves display the same general features as the data for $`G_{Ep}`$, $`G_{Mp}`$ and $`G_{Mn}`$: for low $`Q^2`$ the form factor is close to the dipole form while for large $`Q^2`$ one finds an asymptotic limit for $`G/G_D`$ that depends sensitively upon $`b`$ but is less than unity for reasonable values. The greatest sensitivity to the shape of the density is found in a transition region for $`Q^2`$ that ranges from several tenths to several (GeV/$`c`$)2, depending upon $`b`$. Thus, data with similar general features can be fit by modulating a basic Gaussian with an even polynomial, where the polynomial degree can be minimized by an optimal choice of $`b`$. For $`G_{En}`$ one need only require the polynomial part of $`\stackrel{~}{\rho }(k)`$ to begin with $`k^2`$ to ensure that the net charge vanishes. Expansions of this form are no more complicated than other parametrizations in common use, but are free of unphysical cusps at origin.

## III Analysis procedures

### A Linear expansions

To extract radial densities from the nucleon form factor data we employ techniques originally developed for fitting radial distributions to data for scattering of electrons or protons from nuclei . Simple models with a small number of parameters do not offer sufficient flexibility to provide a realistic estimate of the uncertainty in a radial density. Rather, we employ linear expansions in complete sets of basis functions that are capable of describing any plausible radial distribution without strong a priori constraints upon its shape. Such expansions permit one to estimate the uncertainties in the fitted density due to both the statistical quality of the data and the inevitable limitation of experimental data to a frequency range, $`kk_{max}`$, where

$$k_{max}=\frac{Q_{max}}{\sqrt{1+\frac{Q_{max}^2}{4m^2}}}$$
(35)
is the maximum spatial frequency sampled by experimental data limited to $`QQ_{max}`$. The uncertainty due to limitation of $`k`$ is known as incompleteness error.

A radial density can be represented as an expansion of the form

$$\rho (r)=\underset{n}{}a_nf_n(r)$$
(36)
where the basis functions $`f_n(r)`$ are drawn from any convenient complete set. The corresponding Fourier transform then takes the form

$$\stackrel{~}{\rho }(k)=\underset{n}{}a_n\stackrel{~}{f}_n(k)$$
(37)
where

$$\stackrel{~}{f}_n(k)=_0^{\mathrm{}}𝑑rr^2j_0(kr)f_n(r)$$
(38)
represent basis functions in momentum space. The expansion coefficients, $`a_n`$, are fitted to form factor data subject to several minimally restrictive constraints to be discussed shortly. Analyses of this type are often described as model independent because a complete basis can reproduce any physically reasonable density; if a sufficient number of terms are included in the fitting procedure the dependence of the fitted density upon the assumptions of the model is minimized. By contrast, simple parametrizations like the Galster model severely constrain the shape of the fitted density.

We consider two bases that have been found useful in the analysis of electron or proton scattering data. The present discussion is limited to monopole densities, but generalizations to higher angular momenta are discussed in Refs. .

The Fourier-Bessel expansion (FBE) employs basis functions of the form

$`f_n(r)`$ $`=`$ $`j_0(k_nr)\mathrm{\Theta }(R_{max}r)`$ (40)
$`\stackrel{~}{f}_n(k)`$ $`=`$ $`{\displaystyle \frac{()^nR_{max}}{k^2k_n^2}}j_0(kR_{max})`$ (41)
where $`\mathrm{\Theta }`$ is the unit step function, $`R_{max}`$ is the expansion radius, and $`k_n=n\pi /R_{max}`$. One advantage of the FBE is that the contribution of each term to the form factor is concentrated around its characteristic frequency $`k_n`$ so that a coefficient $`a_n`$ is largely determined by data with $`kk_n`$. The larger the expansion radius $`R_{max}`$, the smaller the spacing between successive $`k_n`$ and the greater the sensitivity one has to variations in the form factor. One should choose $`R_{max}`$ to be several times the root-mean-square radius but not so large that an excessive number of terms is needed to span the experimental range of momentum transfer. Terms with $`k_n>k_{max}`$ provide an estimate of the incompleteness error. We chose $`R_{max}=4.0`$ fm, but the results are insensitive to its exact value. However, a disadvantage of the FBE is that a relatively large number of terms is often needed to accurately represent a typical confined density.

Alternatively, the Laguerre-Gaussian expansion (LGE) employs basis functions of the form

$`f_n(r)`$ $`=`$ $`e^{x^2}L_n^{1/2}(2x^2)`$ (43)
$`\stackrel{~}{f}_n(k)`$ $`=`$ $`{\displaystyle \frac{\sqrt{\pi }}{4}}b^3()^ne^{y^2}L_n^{1/2}(2y^2)`$ (44)
where $`x=r/b`$, $`y=kb/2`$, and $`L_n^a`$ is a generalized Laguerre polynomial. A significant advantage of the LGE is that the number of terms needed to provide a reasonable approximation to the density can be minimized by choosing $`b`$ in accordance with the natural radial scale. We chose $`b=0.556`$ fm such that the mean-square radius of the Gaussian factor is consistent with that of the common dipole parametrization of Sachs form factors. We then find that the magnitude of $`a_n`$ decreases rapidly with $`n`$, but the quality of the fit and the shape of the density are actually independent of $`b`$ over a wide range. However, a disadvantage of the LGE is that the basis functions are not localized in momentum space so that the coefficients tend to be correlated more strongly than for the FBE.

### B Constraints

The expansion coefficients are obtained by minimizing

$$\chi ^2=\underset{i}{}\left(\frac{y_i\overline{y}_i}{\delta y_i}\right)^2$$
(45)
where $`\overline{y}_i`$ is the fitted value of a quantity $`y_i`$ with uncertainty $`\delta y_i`$. In addition to experimental data, the set $`y_i`$ generally includes pseudodata used to enforce constraints and to estimate the incompleteness error associated with the limitation of experimental data to a finite range of momentum transfer.

The absence of data for very large $`Q^2`$ requires some constraint upon the behavior of $`\stackrel{~}{\rho }(k)`$ for $`k_{max}<k<k_m`$. Furthermore, inversion of the Fourier transform also requires an assumption about the experimentally inaccessible region $`k>k_m`$. On quite general grounds one expects the asymptotic form factor for a confined system to decrease more rapidly than $`k^4`$ ; in particular, this condition ensures that there will be no cusp at the orgin. In fact, our results show that the intrinsic form factor $`\stackrel{~}{\rho }(k)`$ is well approximated by a Gaussian for large $`k`$. Therefore, we will assume that $`\stackrel{~}{\rho }(k)`$ for $`k>k_m`$ is bounded by a $`k^4`$ envelope and use the flexibility afforded by that envelope to estimate the incompleteness error due to the limitation of experimental information to the range $`k<k_{max}`$. Although some restriction is needed to stabilize the fits, the $`k^4`$ envelope probably overestimates the uncertainties in unmeasured form factors and their effect upon uncertainties in fitted densities; nevertheless, we prefer to employ minimally restrictive constraints so that those densities will have the best possible model independence.

More detailed discussions of the method may be found in Refs. , but the basic idea is to supplement the experimental data by pseudodata of the form $`\stackrel{~}{\rho }(k_i)=0\pm \delta \stackrel{~}{\rho }(k_i)`$ whose uncertainties are based upon a reasonable model of the asymptotic behavior of the form factor for $`k_i>k_{max}`$ where $`k_{max}`$ is the spatial frequency corresponding to the maximum measured $`Q^2`$. Therefore, uncertainties in the form factor for $`k>k_{max}`$ are based upon an envelope of the form

$`\delta \stackrel{~}{\rho }(k)`$ $`=`$ $`\sqrt{{\displaystyle \frac{1}{3}}}\rho _{lim}(k)`$ (47)
$`\rho _{lim}(k)`$ $`=`$ $`\left|\stackrel{~}{\rho }(k_{max})\right|\left({\displaystyle \frac{k_{max}}{k}}\right)^4`$ (48)
where the factor of $`\sqrt{1/3}`$ represents the variance of a uniform distribution of unit width. When using FBE the pseudodata are chosen at the characteristic frequencies $`k_n=n\pi /R_{max}`$ with $`n>k_{max}R_{max}/\pi `$, while a uniform spacing of $`\mathrm{\Delta }k_i=1.0`$ fm-1 was employed for LGE. The error band for a fitted density is computed from the covariance matrix for the $`\chi ^2`$ fit and includes the incompleteness error. A detailed discussion of the decomposition of the density uncertainty into statistical and incompleteness errors may be found in Ref. .

Recognizing that pQCD imposes an asymptotic limit of the form $`GQ^4`$ upon the Sachs form factors, one might be tempted to employ pseudodata for $`G`$ at large $`Q^2`$. If one knew how to estimate the proportionality constant, this procedure could be used to regulate $`\stackrel{~}{\rho }(k)`$ for $`k_{max}<k<k_m`$ but would not be sufficient for construction of the radial density because inversion of the Fourier transform also requires information for the $`k>k_m`$ region that is inaccessible to electron scattering. Although we expect $`\stackrel{~}{\rho }(k)`$ to be small for $`k>2m`$, we cannot simply set it to zero because an abrupt cut-off would introduce unreasonable density oscillations at very large radii. The present procedure estimates the uncertainty in the radial density arising from both the unmeasured and the unmeasurable ranges of spatial frequency. In this model, the minimum uncertainty in density is governed by the nucleon Compton wavelength and can be interpreted as an irreducible smearing by zitterbewegung.

Small but undesirable oscillations in fitted densities at large radii were suppressed using a tail bias based upon the method discussed in Ref. . We employed a tail function of the form $`t(r)e^{\mathrm{\Lambda }r}`$, based upon the successful dipole parametrization for low $`Q^2`$, and included in the $`\chi ^2`$ fit a penalty function of the form

$$\chi _r^2=\underset{i=1}{\overset{N_d}{}}\left(\frac{\rho (r_i)t(r_i)}{wt(r_i)}\right)^2$$
(49)
to suppress strong deviations from the tail function. The radial pseudodata were constucted on the grid $`r_i=r_m+i\mathrm{\Delta }r`$ for $`i=1,N_d`$ in the range $`r>r_m`$. We chose $`\mathrm{\Lambda }=4.27`$ fm-1, $`r_m=2.0`$ fm, $`\mathrm{\Delta }r=0.2`$ fm, $`N_d=10`$, and $`w=2`$, but the results are rather insensitive to these details. The tail bias improves the convergence of moments of the density but has little effect upon a fitted density in the region where it is large.

The fitting procedure also permits constraints to be placed on integral moments of the radial density. We define fitted moments by

$$\overline{M}_\alpha =_0^{\mathrm{}}𝑑rr^{2+\alpha }\rho (r)$$
(50)
where here $`\rho (r)`$ is the fitted density and include a penalty function of the form

$$\chi _M^2=\left(\frac{M_0\overline{M}_0}{\delta M_0}\right)^2+\left(\frac{M_2\overline{M}_2}{\delta M_2}\right)^2$$
(51)
where $`M_\alpha `$ is the measured value and $`\delta M_\alpha `$ is its uncertainty. The constraint on the neutron charge was enforced by means of a pseudodatum $`M_0=0\pm 10^6`$. In addition, the atomic physics datum for $`M_2`$ from Ref.  was included in fits made to the neutron charge density.

It is also useful to define a fitted transition radius $`\overline{R}`$ as

$$\overline{R}=\sqrt{\overline{M}_2/\overline{M}_0}$$
(52)
for $`G_{Ep}`$, $`G_{Mp}`$, or $`G_{Mn}`$. Thus, the fitted transition radius is correlated with the experimental normalization at low $`Q^2`$. If the fitted $`\overline{M}_0`$ were constrained, the uncertainty in $`\overline{R}`$ would be artificially reduced. Therefore, no constraints were placed on $`M_0`$ for $`G_{Ep}`$, $`G_{Mp}`$, or $`G_{Mn}`$ — given that those intrinsic densities were defined with unit normalization, the fitted values of $`\overline{M}_0`$ test the normalization of the experimental data.

### C Data selection

We tried to select the best available data in each range of $`Q^2`$, with an emphasis upon recent data using recoil or target polarization wherever available. These selections are summarized in Table I. Although a thorough review of the data for nucleon electromagnetic form factors is beyond the scope of the present work, in this section we provide brief explanations for some of our selections and omissions.

$`G_{Mp}`$ data were taken from the compilation of Höhler et al for $`Q^2<0.15`$ (GeV/$`c`$)2 and for larger $`Q^2`$ from the compilation and re-analysis made by Brash et al to improve the correction for the small contribution of $`G_{Ep}`$ to cross section at large $`Q^2`$. Values for $`G_{Ep}`$ were obtained by multiplying the recent recoil polarization measurements of $`G_{Ep}/G_{Mp}`$ from Refs.  by the Brash parametrization of $`G_{Mp}`$. Supplementary recoil polarization data from Ref.  were omitted — those data are consistent with those selected but have larger statistical uncertainties. In addition, the cross section data for $`G_{Ep}`$ from Refs.  were used at low $`Q^2`$. The cross section data for $`G_{Ep}`$ at $`Q^2>1`$ (GeV/$`c`$)2 from Refs.  were omitted because, as shown in Ref. , they are in significant disagreement both with the recoil polarization data and with each other presumably because the Rosenbluth technique becomes increasingly difficult as $`Q^2`$ increases and the relative contribution of the electric form factor to the unpolarized cross section becomes quite small. Nevertheless, we eagerly await new results from a proposed improvement of the Rosenbluth method .

The neutron magnetic contribution is large enough at high $`Q^2`$ to employ quasielastic electron-deuteron scattering with subtraction of the proton contribution. Data of this type were obtained from Refs. . At low $`Q^2`$ the model dependence of the quasielastic method becomes relatively large. Markowitz et al.  measured the quasifree neutron knockout cross section for the $`d(e,e^{}n)`$ reaction and calibrated the efficiency of the neutron detector using associated particle production in the deuteron photodisintegration reaction, $`d(\gamma ,pn)`$. The dependence upon the deuteron wave function can be reduced by analyzing the ratio between quasifree cross sections for neutron or proton knockout, indicated as the $`d(e,e^{}N)`$ reaction in Table I, with relatively small corrections made for meson-exchange currents and final-state interactions (FSI). Bruins et al.  calibrated their neutron detector using the $`p(\gamma ,\pi ^+n)`$ reaction while Refs.  employed elastic neutron-proton scattering. The associated-particle techniques permit calibration in situ but must correct the bremsstrahlung measurements for the contributions of three-body reactions from electroproduction that lie outside the acceptance , whereas the $`p(n,p)n`$ reaction is kinematically complete but requires calibration at a different facility and under different conditions than used for the reaction of interest. Unfortunately, these methods remain in substantial disagreement; because we are not convinced there is a compelling preference, we include the data from both methods in the present analysis. Finally, we also include recent low $`Q^2`$ data  from inclusive electron scattering from transversely polarized $`{}_{}{}^{3}\stackrel{}{\mathrm{He}}`$ that uses the Fadeev calculations of Refs.  to correct for nuclear structure effects. These data are more consistent with the coincidence ratio method calibrated by elastic scattering than by associated-particle production.

We do not use any elastic or quasielastic cross section data for $`G_{En}`$ because the uncertainties arising from nuclear structure are prohibitively large. Polarization techniques offer a signal that is linear in $`G_{En}`$ and with less model dependence. Nevetheless, at low $`Q^2`$ it remains important to correct recoil polarization data for final-state interactions and target polarization for nuclear structure. Most of the data for deuterium targets have been analyzed using the calculations of Arenhövel et al.  to correct for nuclear structure; in Table I we cite both the experimental paper and the subsequent analysis. The result from the first experiment of this type  has not been corrected, but the statistical uncertainty was large. The $`Q^2=0.4`$ (GeV/$`c`$)2 data for $`{}_{}{}^{3}\stackrel{}{\mathrm{He}}(\stackrel{}{e},e^{}n)`$ have been analyzed using Fadeev calculations  but Rohe et al.  argue that the plane-wave impulse approximation (PWIA) is adequate at larger $`Q^2`$. We also include values for $`G_{En}`$ extracted by Schiavilla and Sick  from an analysis of the deutron quadrupole form factor obtained from tensor polarization measurements of elastic scattering . Although the model uncertainties remain somewhat large, this analysis covers a larger range of $`Q^2`$ and appears to be more accurate than the older analysis of deuteron elastic scattering by Platchkov et al. , which is omitted.

Finally, the neutron mean-square charge radius is related to the neutron-electron scattering length, $`b_{ne}`$, by

$$r^2_n=\frac{3\mathrm{}}{\alpha m_nc}b_{ne}$$
(53)
Unfortunately, the measurements are rather difficult and most techniques require substantial corrections for effects for which there is often insufficient information. A recent review of these measurements has been made by Alexandrov , who finds that most modern measurements cluster around two values. From measurements of the energy dependence for the transmission of thermal neutrons through liquid 208Pb, Kopecky et al. Ref.  obtained $`b_{ne}=(1.33\pm 0.027\pm 0.03)\times 10^3`$ fm, corresponding to $`r^2_n=0.115\pm 0.003`$ fm2. This result agrees well with similar measurements by Koester et al.  for lead isotopes and 209Bi and with the results of Krohn and Ringo  using the angular distribution for neutron scattering by noble gases. Alternatively, Alexandrov et al.  obtained $`b_{ne}=(1.60\pm 0.05)\times 10^3`$ fm, corresponding to $`r^2_n=0.138\pm 0.004`$ fm2, using neutron diffraction from single crystals of 186W. This result is consistent with a bismuth transmission experiment that was also performed at Dubna, but disagrees by about 5 standard deviations from the Garching, Argonne, and Oak Ridge experiments. This discrepancy has been attributed to resonance corrections  but remains controversial. The extracted scattering length is strongly correlated with the resonance correction. Leeb and Teichtmeister  argue that the correction employed by Alexandrov et al.  requires implausibly large contributions from negative energy levels. Alexandrov  argues that the energy-independent resonance correction should be fitted to the data and that in the absence of definitive information a negative contribution cannot be excluded. Futhermore, neutron diffraction from single crystals of 186W provides a larger signal than the energy dependence of the total cross section, and is hence is less sensitive to this correction. Although this controversy has not yet been resolved satisfactorily, we decided to employ the most recent result from Oak Ridge, namely $`r^2_n=0.115\pm 0.003`$ fm2, as a datum in our fit of the neutron charge density and to omit the Dubna result. The sensitivity to this choice is discussed in Sec. IV C 4.

## IV Results

### A Form Factors

Fits to the form factor data are shown in Fig. 2 as bands that represent the uncertainties in the fitted form factors. These bands were computed using the covariance matrix. The fits shown in Fig. 2 employ the LGE parametrization with $`\lambda _E=\lambda _M=2`$, but the results using the FBE parametrization are practically indistinguishable. Nor do these fits depend upon the choices for $`b`$, $`R_{max}`$, or details of the constraints. Fits using $`\lambda <2`$ are almost identical within the ranges spanned by experimental data, but their error bands grow more rapidly at larger $`Q^2`$. The rapidly decreasing dipole form factor is divided out to emphasize the deviations at large $`Q^2`$ from this characteristic behavior. For $`G_{En}`$ we also display a simple two-parameter fit using the Galster parametrization.

The intrinsic form factors, obtained via the relativistic transformation prescribed by Eq. (II B), are shown in Fig. 3 using $`\lambda _E=\lambda _M=2`$. From these figures we observe that for moderate $`k^2`$ three of the four intrinsic form factors resemble simple Gaussians, while the intrinsic neutron charge form factor requires an additional factor of $`k^2`$ in first approximation. Consequently, only a few terms of the Laguerre-Gaussian expansion are needed to obtain good fits, with higher-order terms used primarily for the estimation of the incompleteness error. Although it is possible to obtain fairly good fits using just 2 terms for $`G_{En}`$, 4 for $`G_{Ep}`$, or 6 for $`G_{Mp}`$ and $`G_{Mn}`$, in order to minimize model dependence and to evaluate incompleteness errors we employed 20 terms for each of the four form factors.

The widths of the form factor bands are governed by the quality of the experimental data in the ranges of $`k^2`$ where data are available and for larger $`k^2`$ by the asymptotic envelopes indicated by dashed curves. Note that the uncertainties in the fitted form factors for $`k>k_{max}`$ are reduced by the factor of $`1/3`$ used in Eq. (III B) to transform from a uniform to a normal distribution and by the effect of constraints upon densities at large radii. Although the intrinsic form factors fitted to data for $`G_{Ep}`$, $`G_{Mp}`$, and $`G_{Mn}`$ appear to decrease more rapidly than the $`k^4`$ envelopes, we prefer to employ these more generous uncertainties rather than to impose the steeper declines suggested by extrapolation from the measured into the unmeasurable region (where $`k>k_m`$). The use of steeper envelopes would simply reduce the uncertainties in the extracted densities without affecting their central values. Therefore, $`k^4`$ envelopes were matched to fitted form factors at $`k_{max}=6.5`$, 9.0, and 8.2 fm-1 for $`G_{Ep}`$, $`G_{Mp}`$, and $`G_{Mn}`$ based upon the experimental $`Q_{max}`$ for each form factor. However, uncritical application of the same procedure to $`G_{En}`$ would suppress the high-frequency components of its intrinsic form factor too strongly because $`Q_{max}`$ for $`G_{En}`$ is presently too small to expect $`\stackrel{~}{\rho }(k)`$ to decrease more rapidly than $`k^4`$. Fig. 2 shows that the data presently available for $`G_{En}`$ are compatible with the Galster parametrization, but the procedure used for the other form factors would cause $`G_{En}/G_D`$ to decrease fairly rapidly beyond the range of these data. On the other hand, it is reasonable to expect $`G_{En}/G_D`$ to decrease for $`Q^2`$ beyond a few (GeV/$`c`$)2, as observed in the other form factors. Therefore, in order to permit the positive slope for $`G_{En}/G_D`$ to continue over a limited but larger range of $`Q^2`$, we used the same value of $`k_{max}`$ for both $`G_{En}`$ and $`G_{Ep}`$ even though the $`G_{En}`$ data are limited to $`k<5.4`$ fm-1. We believe that this compromise provides a more reliable extrapolation to higher $`Q^2`$ and that the increased estimate of incompleteness error is more realistic, but obviously it is very important to acquire accurate data for $`G_{En}`$ at higher $`Q^2`$.

These fits to intrinsic magnetic form factors do not change sign within the experimentally accessible region, $`k<k_m`$, but the fitted proton intrinsic charge form factor suggests a node between the present experimental limit, $`k_{max}=6.5`$ fm-1, and ultimate limit, $`k_m=9.5`$ fm-1. Consequently, this model suggests a zero crossing in $`G_{Ep}`$ near $`Q^210`$ (GeV/$`c`$)2. Figure 4 compares the form factor ratio $`\mu G_E/G_M`$ deduced from the fitted form factors with experimental data using recoil polarization for the proton or using either recoil or target polarization for the neutron. For the proton we also show the linear parametrization proposed by Jones et al.  for $`Q^2>0.3`$ (GeV/$`c`$)2, while for the neutron we show a new fit using the Galster parametrization, Eq. (2), that gave $`A=0.90\pm 0.02`$ and $`B=3.8\pm 0.5`$. The data for the proton do not distinguish between linear and LGE parametrizations, but according to pQCD one would expect $`G_{Ep}/G_{Mp}`$ to approach a constant for sufficiently large $`Q^2`$. Extrapolation of the LGE parametrization suggests that the asymptotic ratio will be very small, but data at much larger $`Q^2`$ are needed to establish that level. An extension to 9 (GeV/$`c`$)2 has been approved , but larger $`Q^2`$ remains desirable. Similarly, the present data for $`G_{En}/G_{Mn}`$ are compatible with the Galster parametrization but remain limited to rather small $`Q^2`$. Consequently, the extrapolation to larger $`Q^2`$ is rather uncertain. If an approved experiment using the $`{}_{}{}^{3}\stackrel{}{\mathrm{He}}(\stackrel{}{e},e^{}n)`$ reaction  achieves the proposed $`\pm 13\%`$ statistical uncertainty at $`Q^2=3.4`$ (GeV/$`c`$)2, the error band will be reduced to about the same width and the extrapolation much improved. Nevertheless, there is little reason to expect the asymptotic limit to be reached earlier for the neutron than for the proton.

Although a review of recent theoretical calculations is beyond the scope of the present work, it is probably worth mentioning a few which describe the new $`G_{Ep}/G_{Mp}`$ data relatively well. Among these the earliest is the chiral soliton model of Holzwarth , which predicted the linear descrease with respect to $`Q^2`$ and a sign change near 10 (GeV/$`c`$)2. More recently , modifications of the vector meson parameters were made to improve the fits to the neutron form factors, but the ratio $`G_{Mn}/G_{Mp}`$ is not reproduced. Furthermore, because the chiral soliton model uses $`\lambda _E=0`$ and $`\lambda _M=1`$, Holzwarth found it necessary to artificially increase the soliton mass in order to obtain reasonable fits at large $`Q^2`$. Alternatively, Lu et al.  obtained a good fit to the $`G_{Ep}/G_{Mp}`$ data for $`Q^23`$ (GeV/$`c`$)2 by adjusting the bag radius in the cloudy bag model, but the ratio appears to level off well above the more recent data for higher $`Q^2`$. Note that this model uses $`\lambda _E=\lambda _M=1`$. The covariant calculation of Boffi et al.  using the point-form spectator approximation provides reasonably accurate predictions of the form factors for $`Q^25`$ (GeV/$`c`$)2, although there remains a significant discrepancy for $`G_{Mp}`$ near the end of this range. The light-front calculations of Cardarelli and Simula  using one-gluon exchange and the light-cone diquark model of Ma et al.  also reproduce the linear $`Q^2`$ dependence of $`G_{Ep}/G_{Mp}`$ fairly well.

### B Densities

Proton charge and magnetization densities are compared in Fig. 5. Both densities are measured very precisely, with uncertainties at the origin better than 6% for magnetization or 8% for charge. Incompleteness dominates in the interior region while statistical errors become comparable in the surface region. As shown by the variation of $`G_{Ep}/G_{Mp}`$ in the top panel of Fig. 4, the new recoil-polarization data for $`G_{Ep}`$ decrease more rapidly than either the dipole form factor or the magnetic form factor for $`Q^2>1`$ (GeV/$`c`$)2. Consequently, we find that the charge density is significantly softer than the magnetization density of the proton. The densities obtained using LGE or FBE parametrizations are practically indistinguishable and are independent of the choice of $`b`$ or $`R_{max}`$ over wide ranges. These densities are similar to the Gaussian densities one might expect from a quark model and are more realistic than the exponential density that results from naive nonrelativistic inversion of the dipole form factor.

Neutron densities are shown in Fig. 6. We find that the magnetization density for the neutron is very similar to that for the proton, although the interior precision is not as good because the range of $`Q^2`$ is smaller and the experimental uncertainties larger. Limitations in the range and quality of the $`G_{En}`$ data presently available result in a substantially wider error band for the neutron charge density. Data at higher $`Q^2`$ are needed to improve the interior precision, but a useful measurement of the interior charge density is obtained nonetheless. The positive interior density is balanced by a negative surface lobe. Note that polarization measurements are sensitive to the sign of the density.

Whereas Figs. 5-6 emphasize the interior densities, it is also of interest to compare these densities in the surface and tail regions. Figures 7-8 use a factor of $`r^2`$ to emphasize these surface and tail densities. Although the densities are small, the reduced slopes seen between 1 and 1.5 fm in the neutron magnetization and in both the charge and the magnetization densities for the proton are seen as significant peaks in $`r^2\rho `$. Virtually identical features also emerge using the FBE parametrization. These features are independent of $`b`$ for the LGE or $`R_{max}`$ for the FBE parametrization over wide ranges. Many attempts were made to suppress structures in $`r^2\rho `$ in the 1-1.5 fm region by limiting the number terms in the expansions or by application of stricter tail biases, but all modifications which did produce smoother $`r^2\rho `$ curves in this region damaged the fits to the form factor data for $`Q^21`$ (GeV/$`c`$)2. Although it is difficult to prove that smoother fits do not exist, especially if one is willing to tolerate a moderate increase in $`\chi ^2`$, we were unable to produce acceptable fits without some structure in $`r^2\rho `$ in this region. On other hand, because local errors in momentum space can introduce artificial oscillations at large radii, we did apply an exponential tail bias for $`r>2`$ fm where little structure is expected. Thus, the smaller oscillations for $`r>1.5`$ fm are generally consistent with zero and can be suppressed using the tail bias with little effect upon the fits. We believe that the 2 fm matching radius is sufficiently large to have minimal influence upon densities for intermediate distances governed by data with $`Q^2`$ of order several (GeV/$`c`$)2.

The relatively small differences between $`G_{Mn}`$ and $`G_{Mp}`$ seen in Fig. 2 produce the small differences between neutron and proton magnetization densities shown in Fig. 8. The peak of $`r^2\rho _m`$ is found at a slightly larger radius for the neutron than for the proton because the form factor decreases a little more rapidly with respect to $`Q^2`$. The secondary peaks in the 1-1.5 fm region are also similar. Again, these comparisons are independent of the details of the analysis and are virtually identical using either LGE or FBE parametrizations.

Therefore, the secondary $`r^2\rho `$ peaks in the 1-1.5 fm region appear to be essential features of the data rather than artifacts of analyses based upon linear expansions. While it is not possible to determine the physical mechanism for such features from data analysis alone, there is at least one simple candidate. The tensor interaction between quarks is expected to produce a small $`D`$-state component peaking at larger radius than the dominant $`S`$-state configuration, and the superposition of these components could yield a secondary peak at relatively large radius.

### C Fitted Moments

Moments and $`\chi ^2`$ for each fit are listed in Table II. The expansion coefficients require too much space to list but are available by request; also note that accurate reproduction of the error bands would require full covariance matrices. Here we quote $`\chi ^2/N`$, where $`N`$ is the number of data points, because counting the number of degrees of freedom is not so clear when both high-$`k`$ and large-$`r`$ constraints are applied. For each form factor we find that the 6 fits obtained for 3 possible choices of $`\lambda `$ used in either LGE or FBE parametrizations are essentially identical and give the same values of $`\chi ^2/N`$. The normalization for $`G_{Ep}`$ is consistent with unity within the 0.5% systematic uncertainty claimed by Simon et al.  for their data at low $`Q^2`$. For the $`G_{Mp}`$ at low $`Q^2`$ we employed the results of Höhler et al.  who adjusted the relative normalizations of several data sets to a common standard. The normalization produced by the present fit is consistent with the systematic uncertainty in that standard. Except for $`G_{Mn}`$, data selected from several sources appear to be mutually consistent and the quality of the fitted form factors is very good. Although the low-$`Q^2`$ data for $`G_{Mn}`$ have improved in recent years, significant systematic discrepancies remain. Recent data from Refs.  with small statistical uncertainties suggest a small dip near 0.2 and a peak near 1 (GeV/$`c`$)2. However, the data from Refs.  are inconsistent with the fit and inflate $`\chi ^2`$. Nor do the data for $`G_{Mn}`$ reach sufficiently low $`Q^2`$ to strongly constrain the normalization and the data sets are not entirely consistent either; consequently, the present analysis suggests that a 2% normalization error remains. We decided to retain the data from Refs. , despite their deviation from the fit because we are not entirely convinced of the transportability of the efficiency calibration from a hadron facility to an electron facility.

The final two columns of Table II list the Sachs radius, $`\xi `$, defined by Eq. (32) in terms of the initial logarithmic derivative of the Sachs form factor and the transition radius, $`\overline{R}`$, obtained from integral moments of the density fitted for a specific $`\lambda `$ according to Eq. (52). The dependencies of rms radii for $`G_{Ep}`$, $`G_{Mp}`$, and $`G_{Mn}`$ upon $`\lambda `$ are consistent with Eq. (33), showing a significant discrete ambiguity arising from the model dependence of the form factor to density inversion. By contrast, for each of these form factors all 6 determinations of $`\xi `$ are consistent with each other, demonstrating that the fitted $`\xi `$ is a model-independent property of the Sachs form factor.

#### 1 Proton charge radius

The analysis of Simon et al.  has been accepted for about two decades as the definitive determination of the proton charge radius, but more recently there has been renewed interest in that quantity now that the finite-size corrections have become the dominant uncertainty in theoretical calculations of the 1S Lamb shift in hydrogen. For example, Melnikov and van Ritbergen  argue that the Lamb shift provides the most accurate measurement of the proton charge radius and deduced a value 0.883(14) fm that is somewhat larger than the 0.862(12) fm obtained by Simon et al. Although our definition for the intrinsic charge radius depends upon the choice of $`\lambda _E`$ employed to fit the Sachs form factor, the definition generally employed by QED theorists corresponds to the quantity we labelled as $`\xi _{Ep}`$ that is based upon the initial slope of $`G_{Ep}(Q^2)`$ and is independent of $`\lambda _E`$. Thus, it is the fitted value of $`\xi _{Ep}`$ that should be compared with the Lamb shift result.

Our fit using the LGE parametrization is compared with the data from Simon et al. in Fig. 9. Also shown is the monopole fit made by Simon et al. to the data for $`Q^2<2.3`$ (GeV/$`c`$)2 that were available at that time. That analysis gave a smaller value, $`\xi _{Ep}=0.862\pm 0.012`$ fm, but does not fit the low $`Q^2`$ data as well as our LGE fit. Our fit employs the entire data set described in Sec. III C even though Fig. 9 shows only the lowest $`Q^2`$ region. Fits made using different values of $`\lambda _E`$ or using the FBE parametrization are indistinguishable and give values for $`\xi _{Ep}`$ that are consistent within their quoted uncertainties. Therefore, we claim that $`\xi _{Ep}=0.88\pm 0.01`$ fm represents a model-independent property of the experimental data even if its interpretation as a charge radius depends upon the choice of $`\lambda _E`$. This value is consistent with the rms radius derived by Melnikov and van Ritbergen  from the Lamb shift.

Coulomb distortion may also affect the charge radius obtained by electron scattering. Rosenfelder  analyzed this effect using a distorted wave Born calculation to correct measured electron scattering cross sections for Coulomb distortion, thereby obtaining effectively plane-wave cross sections. The charge radius was then obtained by fitting the resulting form factors with low-order polynomials. Although the form factor corrections were typically less than 1%, the value of $`\xi _{Ep}`$ extracted from the adjusted data of Simon et al. nevertheless increases by about 0.008 – 0.013 fm depending upon the fitting strategy and the degree of the polynomial. However, that analysis employed a Coulomb potential obtained from a charge density of the form

$$\rho (r)=\frac{d^3k}{(2\pi )^3}e^{i𝐤𝐫}\frac{G_E(k^2)}{\sqrt{1+\frac{k^2}{4m^2}}}$$
(54)
This form not only lacks the Lorentz-contraction factor for spatial frequency, but uses a value for $`\lambda _E=1/2`$ that is inconsistent with the high $`Q^2`$ behavior of the form factor. The question of Coulomb distortion of the low $`Q^2`$ form factors merits further investigation, but a consistent relativistic relationship between form factor and density must be employed. Nevertheless, the magnitude of that correction appears to be smaller than the uncertainty in the fitted quantity.

#### 2 Proton magnetization radius

With $`\xi _{Mp}=0.85\pm 0.03`$ fm the proton magnetization radius appears to be slightly smaller than its charge radius, as expected from Fig. 5, but the uncertainty is as large as the difference because the data at very low $`Q^2`$ are less precise for $`G_{Mp}`$ than for $`G_{Ep}`$. A substantial part of the uncertainty in $`\xi _{Mp}`$ is due to the uncertainty in normalization. If we constrain $`M_0`$ to unity, the changes in fitted form factors and densities are relatively small, but we obtain values for $`\xi _{Ep}=0.862\pm 0.006`$ fm and $`\xi _{Mp}=0.835\pm 0.006`$ fm that appear to be much more precise. The constrained fit to $`G_{Ep}`$ is close to the result of Simon et al. shown in Fig. 9 at very low $`Q^2`$. This analysis demonstrates that there is an appreciable difference between the proton charge and magnetization densities, but it also highlights the importance of precise absolute normalization at very low $`Q^2`$.

#### 3 Neutron magnetization radius

Our fits to the $`G_{Mn}`$ data give a value for $`\xi _{Mn}=0.92\pm 0.07`$ fm that has substantial uncertainty because the lack of data for very low $`Q^2`$ permits significant widening of the error band as $`Q^20`$. By contrast, Kubon et al.  obtained a value $`\xi _{Mn}=0.873\pm 0.011`$ fm that appears to be much more precise. However, their continued-fraction parametrization automatically constrains the normalization at $`Q^2=0`$. If we constrain the normalization by requiring $`M_0=1`$, then the LGE analysis gives a value for $`\xi _{Mn}=0.881\pm 0.018`$ fm that agrees with Kubon et al. Note that we included data from Refs.  that were omitted by Kubon et al. and that deviate strongly from our LGE fit, but these data appear to have little influence upon the fitted normalization and radius. Although the fitted normalization is consistent with unity, the 3% uncertainty does have an appreciable effect upon the uncertainty in $`\xi _{Mn}`$. We chose not to constrain the normalization in the standard analysis because the systematic errors in neutron efficiency have been a big problem for $`G_{Mn}`$ measurements and there remains significant scatter among recent experiments at low $`Q^2`$.

#### 4 Neutron charge radius

The charge radius for the neutron can be expressed in terms of Dirac and Pauli form factors as

$$r^2_n=6\frac{dF_{1n}(Q^2)}{dQ^2}|_{Q^20}+\frac{3\kappa _n}{2m^2}$$
(55)
where the first term is sometimes described as the intrinsic radius (e.g. ) while the second term is called the Foldy term and is attributed to a charge separation induced by the zitterbewegung motion of the magnetization density. We prefer to describe the first term as the Dirac radius because it is derived from the initial $`Q^2`$ dependence of the Dirac form factor and to reserve the term intrinsic radius for Eq. (30), which is based upon a moment of the radial function that we identified as the intrinsic charge density. However, for the present purposes it will be clearer to refer to the radius based upon the Sachs form factor as the Sachs radius. The observation that the Foldy term, equal to -0.126 fm2, is by itself almost equal to the mean-square charge radius obtained from $`b_{ne}`$ has generated considerable discussion. Furthermore, as discussed in Sec. III C, a substantial disagreement remains between the results from Dubna and those from Oak Ridge, Garching, and Argonne. For example, Alexandrov  argues that the Foldy term should be discarded and that the pionic cloud should make the mean-square Dirac radius negative. He further claims that experiments giving $`b_{ne}1.31`$ fm are likely to suffer from serious experimental or interpretative errors because the corresponding mean-square Dirac radius would be positive. On the other hand, Glozman and Riska  calculated that the pion loop contribution to the Dirac radius is negligible. Using a Foldy-Wouthuysen analysis of the interaction of a particle with internal structure with an external electromagnetic field, Bawin and Coon  demonstrated that the Foldy term is cancelled by a higher-order term arising from the Dirac form factor, leaving the Sachs radius as the dominant coefficient, independent of the dynamics responsible for the neutron form factors.

From a more microscopic point of view, Isgur used a quark model to argue that the observation of a very small Dirac radius is a potentially misleading accident and that the relativistic boost cancels the Foldy term such that the slope of $`G_{En}`$ does provide the second moment of the intrinsic charge distribution . Cardarelli and Simula  showed that this cancellation depends upon neglect of transverse momenta and that quark spin-spin interactions that break SU(6) symmetry provide a mixed-symmetry $`S^{}`$ component that enhances $`G_{En}`$ and provides a fairly accurate fit to the recoil-polarization data for $`\mu _pG_{Ep}/G_{Mp}`$. Leinweber et al argued that chiral perturbation theory provides model-independent constraints on the dependencies of nucleon magnetic moments and charge radii upon quark masses that demonstrate that the similiarity between the Sachs radius and the Foldy term is purely accidental.

Our fitted mean-square radius for the neutron charge density is largely determined by and is completely consistent with the datum of Ref.  that was included in the analysis. This quantity is free of discrete ambiguity because $`M_0`$ vanishes and the analysis procedure enforced that constraint explicitly. Unfortunately, the electron-scattering data are not sufficiently precise at very low $`Q^2`$ to resolve the controversy concerning the sign of the Dirac radius. In the absence of a datum for $`M_2`$, the fit to the data for $`G_{En}`$ gives $`M_2=0.187\pm 0.04`$ fm, which is about one or two standard deviations larger in absolute magnitude than the Dubna or Oak Ridge results, respectively, but is much less precise than either — the error bands on $`G_{En}`$ with or without the $`M_2`$ datum overlap almost completely, even at low $`Q^2`$. The effect of the constraint on $`M_2`$ is shown in Fig. 10, where the shaded band for $`F_{1n}`$ represents the unconstrained fit to $`G_{En}`$ data that has a negative Dirac radius while the cross-hatched band with a positive Dirac radius includes the $`M_2`$ datum and is consequently narrower as $`Q^20`$. The rather small difference between these fits is confined to $`Q^20.2`$ (GeV/$`c`$)2, but because $`F_{1n}<0`$ over most of its measured range the fit with negative Dirac radius requires a sign change near $`Q^20.14`$ (GeV/$`c`$)2. Therefore, the present electron-scattering data offer very little sensitivity to $`M_2`$. Furthermore, because the nuclear physics corrections needed to extract $`G_{En}`$ for low $`Q^2`$ are substantial even for polarization methods, it is unclear whether one can ever expect better accuracy from nuclear physics than atomic physics measurements of $`M_2`$. We consider neither the theoretical argument for negative Dirac radius nor the limited experimental evidence for a sign change in $`F_{1n}`$ at low $`Q^2`$ compelling.

## V Discussion

### A Comparison with Nonrelativistic Analyses

Several analyses have appeared recently in which the parameters of the Galster model were fitted to selected data for $`G_{En}`$ at low $`Q^2`$ and a density extracted using the nonrelativistic inversion formula given by Eq. (10). Examples of this type are shown in Figs. 11-12. Figure 11 compares fits of this type to recent data using the d$`(\stackrel{}{e},e^{}\stackrel{}{n})`$, $`\stackrel{}{d}(\stackrel{}{e},e^{}n)`$, or $`{}_{}{}^{3}\stackrel{}{\mathrm{He}}(\stackrel{}{e},e^{}n)`$ reactions. The slope of the form factor obtained from the Oak Ridge value for $`b_{ne}`$ is shown as a line segment. The result obtained by Platchkov et al.  from an analysis of elastic cross sections for electron scattering from deuterium using the Paris potential is shown as the dashed curve and lies well below the data obtained from polarization measurements. However, variations of $`\pm 50\%`$ in $`G_{En}`$ were found using different realistic nucleon-nucleon potentials; the result using the Paris potential is quoted most often, but the result using the Argonne V14 potential is closer to the modern data. The dash-dotted curve shows a fit by Schmieden  to the data from Mainz, excluding the point for $`{}_{}{}^{3}\stackrel{}{\mathrm{He}}(\stackrel{}{e},e^{}n)`$ at $`Q^2=0.4`$ (GeV/$`c`$)2; this fit gives the highest values for $`G_{En}`$. The original Galster model is shown as a solid curve and our fit based upon Eq. (2) is shown as a dotted curve and lies between the Paris and Mainz results. The fit by Zhu et al.  is also close to the dotted curve. Note that our fit uses the entire data set described in Sec. III C, including the slope at the origin, while the Mainz fit used a smaller subset. Even though the data employed by Galster et al.  had much larger uncertainties, their result is remarkably close to the present analysis. However, the apparent agreement of the original Galster parametrization with more modern data must be judged as fortuitous.

Buchmann et al argue that the neutron charge density reflects differences between the spatial distributions of constituent quarks induced by the color hyperfine interaction and that because the same interaction is responsible for the $`N\mathrm{\Delta }`$ mass splitting the $`N\mathrm{\Delta }`$ quadrupole form factor is related in a simple manner to $`G_{En}`$. In Ref.  they used the Galster parametrization to fit a selection of the $`G_{En}`$ data and used nonrelativistic inversion to obtain a density similar to the solid curve in Fig. 12.

The nonrelativistic densities obtained from these fits are compared in Fig. 12 to the present relativistic LGE result. Here we chose $`\lambda _E=0`$ to minimize the differences between relativisitic and nonrelativistic inversion formulas; however, it is important to remember that the Galster form factor is inconsistent with the relativistic inversion procedure for $`\lambda _E=0`$. The upper panel, emphasizing the interior density, demonstrates that the naive Fourier transform tends to produce a rather hard core and an unphysical cusp at the origin. The relativistic transformation, by contrast, softens the interior density and eliminates its cusp, providing a much more plausible charge density. The lower panel, emphasizing the surface lobe, shows that all three nonrelativistic densities have positive peaks at smaller radii and less surface charge than our result.

### B Approach to scaling

The asymptotic behavior of the fitted form factors is illustrated in Fig. 13 by multiplication by $`Q^4`$. The uncertainties in the fitted form factors are clearly dominated by experimental uncertainties where data are available, while the expansion of the error bands for larger $`Q^2`$ is governed by the incompleteness errors in the range $`k_{max}<k<k_m`$. The scaling behavior of $`G_{Mp}`$ appears fully developed because the data up to $`Q_{max}^2=31`$ (GeV/$`c`$)2 reach a sufficiently large ratio $`k_{max}/k_m=0.95`$ to strongly constrain the asymptotic limit. Although the fit to $`G_{Mn}`$ is compatible with pQCD scaling, the uncertainties at large $`Q^2`$ grow much more rapidly because here $`k_{max}/k_m=0.86`$ is smaller and the data for $`k_{max}/k_m>0.79`$ are much less precise. Furthermore, because scaling is not fully developed for $`G_{Mp}`$ until $`Q^220`$ (GeV/$`c`$)2 we should expect that data for $`G_{Mn}`$ at higher $`Q^2`$ will be needed to establish its asymptotic limit. The data for $`G_{Ep}`$ do not show scaling behavior for $`Q_{max}^2=5.54`$ (GeV/$`c`$)2, where $`k_{max}/k_m=0.78`$ — despite relatively large uncertainties in extrapolation, the present fit suggests a sign change near $`Q^210`$ (GeV/$`c`$)2. The new recoil polarization data suggest that $`G_{Ep}`$ and $`G_{Mp}`$ differ dramatically for $`Q^23`$ (GeV/$`c`$)2 and it is clearly crucial to extend the $`G_{Ep}`$ data to higher $`Q^2`$. It will also be important to check those data using another technique. The results of an improved Rosenbluth experiment are expected soon .

The present data for $`G_{En}`$ are too limited in both range and precision to address the question of scaling. High-precision data for $`Q^21.5`$ (GeV/$`c`$)2 from a recent d$`(\stackrel{}{e},e^{}\stackrel{}{n})`$ experiment are expected soon , and an approved proposal for $`{}_{}{}^{3}\stackrel{}{\mathrm{He}}(\stackrel{}{e},e^{}n)`$ should extend the range to 3.4 (GeV/$`c`$)2 in a couple of years . However, experience with the other three form factors suggests that one must approach 20 (GeV/$`c`$)2 to determine the asymptotic limit.

Quark helicity conservation suggests that $`Q^2F_2/F_1`$ should approach a constant in the asymptotic limit of large $`Q^2`$. The Dirac and Pauli form factors are related to Sachs form factors by

$`F_1`$ $`=`$ $`{\displaystyle \frac{G_E+\tau G_M}{1+\tau }}`$ (57)
$`\kappa F_2`$ $`=`$ $`{\displaystyle \frac{G_MG_E}{1+\tau }}`$ (58)
and their ratio is given by

$$\frac{\kappa F_2}{F_1}=\frac{1g}{\tau +g}$$
(59)
where $`g=G_E/G_M`$ can be measured directly by either recoil or target polarization. Recoil-polarization data for this ratio are compared in Fig. 14 with bands constructed from the present fits to Sachs form factors. Gayou et al.  observed that the proton recoil polarization data appear to reach a plateau in the range $`2<Q^2<6`$ (GeV/$`c`$)2 when scaled by $`Q`$ instead of the expected $`Q^2`$. Ralston et al.  suggested that orbital angular momentum in the quark distribution could explain an asymptotic behavior of the form $`F_2/F_1Q^1`$. Later Miller and Frank  argued that substantial violation of quark helicity conservation should be expected for intermediate $`Q^2`$ when Poincaré invariance is imposed upon relativistic constituent quark models.

Recognizing that $`g`$ is small compared with $`\tau `$ for $`Q^2>6`$ (GeV/$`c`$)2 and that the model imposes constraints upon the Sachs form factors for high $`Q^2`$ that inhibit the growth of the error band for $`F_2/F_1`$, we can extrapolate $`F_2/F_1`$ beyond the present experimental range for $`G_{Ep}`$. This extrapolation, shown in Fig. 15, shows that the data are consistent with quark helicity conservation for $`Q^220`$ (GeV/$`c`$)2. Therefore, although the present model is consistent with the observation by Gayou et al. that $`QF_2/F_1`$ is approximately constant in the range $`2<Q^2<6`$ (GeV/$`c`$)2, we attribute that observation to the appearance of a broad maximum in $`QF_2/F_1`$ rather than to the onset of true $`Q^1`$ scaling.

### C Quark densities

Assuming isospin symmetry and neglecting strange quarks, nucleon charge densities can be expressed in a simple two-flavor quark model as

$`\rho _p(r)`$ $`=`$ $`{\displaystyle \frac{4}{3}}u(r){\displaystyle \frac{1}{3}}d(r)`$ (61)
$`\rho _n(r)`$ $`=`$ $`{\displaystyle \frac{2}{3}}u(r)+{\displaystyle \frac{2}{3}}d(r)`$ (62)
where $`u(r)`$ is the radial distribution for an up quark in the proton or a down quark in the neutron while $`d(r)`$ is the distribution for a down quark in the proton or an up quark in the neutron. Thus, the quark densities are obtained from nucleon charge densities using

$`u(r)`$ $`=`$ $`\rho _p(r)+{\displaystyle \frac{1}{2}}\rho _n(r)`$ (64)
$`d(r)`$ $`=`$ $`\rho _p(r)+2\rho _n(r)`$ (65)
where $`u(r)`$ and $`d(r)`$ are normalized to unity according to

$$_0^{\mathrm{}}𝑑rr^2q(r)=1$$
(66)
but need not be positive everywhere. There is no guarantee that these combinations of radial densities obtained from form factor data by relativistic inversion must be positive, nor are the densities derived from positive-definite matrix elements. Within the quark model one could decompose the densities

$$q(r)=q_v(r)+q_s(r)\overline{q}_s(r)$$
(67)
for each flavor, $`q`$, into valence ($`v`$) and sea ($`s`$) contributions which would be expected to be positive definite but cannot be separated using only Sachs form factors. The valence distributions, $`q_v`$, are normalized to unity while the sea distributions, $`q_s`$ and $`\overline{q}_s`$, must have equal normalizations. Thus, within this picture the contributions from quarks should be positive but antiquarks in the sea could produce regions where $`u(r)`$ or $`d(r)`$ might become negative.

The quark densities obtained with LGE densities for $`\lambda _E=2`$ that are displayed in Fig. 16 are determined with relatively small uncertainties and are predominantly positive, as expected. We find that the pair of like quarks has a somewhat broader distribution than that of the unlike quark, being depleted in the interior and enhanced at the surface. In this model, the neutron charge density arises from incomplete cancellation between charge densities for up and down quarks, resulting in positive core and negative surface charges. The broader distribution for like quarks is consistent with the repulsive color hyperfine interaction between like quarks needed to explain the $`N\mathrm{\Delta }`$ mass difference. This picture is also consistent with the model of a pion cloud surrounding a three-quark core, but in that model one might expect to find a slightly negative $`d(r)`$ near the surface due to the antiquark content of the pion. The present data are marginally consistent with a slightly negative $`d(r)`$ near 1.0 fm, but more accurate data for the neutron charge density would be needed to reduce the uncertainty in $`d(r)`$ before drawing a definitive conclusion.

### D Discrete Ambiguities

Although discrete ambiguities in $`\lambda _E`$ and $`\lambda _M`$ do not affect fitted form factors in the range where data are available, the choice of $`\lambda `$ does affect the growth of the uncertainties in extrapolated form factors as $`Q^2`$ increases beyond the measured range. For example, Fig. 17 shows form factors fitted using the LGE parametrization with $`\lambda _E=0`$ and $`\lambda _M=1`$ as suggested by the relativistic soliton model. Although the relative uncertainties become quite large for $`Q^2`$ beyond the range of the experimental data, the uncertainties in the form factors actually remain small because, with the exception of $`G_{En}`$, the form factors are greatly reduced as their experimental limits are reached. Consequently, the contribution of the uncertainties in form factors at large $`Q^2`$ upon uncertainties in radial densities is relatively insensitive to the choice of $`\lambda `$. However, the change in fitted density due to a change of $`\lambda `$ need not be contained within the fitted error band — the bands do not accommodate discrete ambiguities in the model.

The sensitivity of fitted charge densities to the choice of $`\lambda _E`$ is illustrated in Fig. 18. These figures were made with the LGE parametrization, but very similar results are obtained with the FBE parameterization. The smoothest results at large radii are obtained with $`\lambda =0`$, whereas larger values of $`\lambda `$ tend to pull the density inward and to amplify oscillations at large radii. This behavior can be understood by interpreting Eq. (II B) in terms of the convolution theorem for Fourier transforms. Expressing this equation in terms of $`k`$,

$$\stackrel{~}{\rho }(k)=(1\frac{k^2}{k_m^2})^\lambda G(\frac{k^2}{1\frac{k^2}{k_m^2}})$$
(68)
one finds that $`\stackrel{~}{\rho }(k)`$ for $`\lambda >0`$ is obtained from the Lorentz contracted form factor by deconvolution of a resolution function with mean square radius equal to $`3\lambda /2m^2`$. This resolution function originates in the zitterbewegung and is characterized by the nucleon Compton wavelength. As discussed in Sec. IV B, acceptable fits to the Sachs form factor data for $`Q^21`$ (GeV/$`c`$)2 using $`\lambda =2`$ seem to require structure in the radial densities in the 1-1.5 fm region. Reducing $`\lambda `$ tends to smooth out such structures, but sacrifices the high $`Q^2`$ limit. Therefore, although $`\lambda =2`$ provides the most natural implementation of pQCD scaling, accurate reproduction of the data using that representation of the intrinsic form factor appears to require oscillations in the radial density with a wavelength of order 0.7 fm.

Despite this ambiguity in the relationship between form factor and density, the qualitative conclusion that $`u(r)`$ is broader than $`d(r)`$ depends only the assumption of isospin symmetry and the observation that the neutron charge density is positive in the interior and negative at its surface. The quark distributions derived using $`\lambda _E=0`$ shown in Fig. 19 are qualitatively similar to those shown for $`\lambda _E=2`$ in Fig. 16, but are slightly more diffuse. The choice $`\lambda _E=0`$ is appropriate for the soliton model, is consistent with nonrelativistic expectations for small $`Q^2`$, and is favored by the radius obtained from the Lamb shift. However, it appears to be inconsistent with pQCD. Therefore, in the absence of a unique relationship between form factors and densities, it appears necessary to select the appropriate value of $`\lambda `$ based upon the intended application. For long-wavelength properties one should use $`\lambda _E=0`$, but for extrapolation to the pQCD limit one should employ $`\lambda _E=\lambda _M=2`$.

### E Alternative parametrizations

Without a unique relationship between form factor and intrinsic density, one may question the value of densities extracted by the present techniques. It is clear that simple parametrizations like the dipole model cannot fit the proton or $`G_{Mn}`$ data over a wide range of $`Q^2`$; nor is it likely that the Galster parametrization will continue to fit $`G_{En}`$ at higher $`Q^2`$. The empirical parametrization proposed by Bosted

$$G(1+a_1Q+a_2Q^2+a_3Q^3+a_4Q^4)^1$$
(69)
fits the data for large $`Q^2`$ well and is consistent with pQCD, but its odd powers of $`Q`$ are incompatible with the interpretation of the form factor as the Fourier transform of a radial density and with the moment expansion for small $`Q^2`$. Furthermore, it is not sufficiently flexible to provide realistic error bands, especially if the odd powers are eliminated. By inclusion of a $`Q^5`$ term, Brash et al.  also sacrificed the pQCD limit in order to improve the quality of the fit for finite $`Q^2`$.

Kubon et al.  fit a subset of the $`G_{Mn}`$ data using a continued-fraction parametrization of the form

$$G_{Mn}(Q^2)=\frac{\mu _n}{1+\frac{b_1Q^2}{1+\frac{b_2Q^2}{1+\mathrm{}}}}$$
(70)
carried to fifth order. This parametrization provides a good fit to the data for $`Q^2<4`$ (GeV/$`c`$)2 using 5 parameters, but the parameters do not decrease with order and the fit depends upon fairly delicate cancellations. Adding additional terms to extend the range of $`Q^2`$ changes the lower terms. Furthermore, this fit does not conform with the asymptotic $`Q^4`$ behavior expected by pQCD unless a fairly complicated constraint of the form

$$b_3b_5+b_2(b_4+b_5)=0$$
(71)
is imposed to eliminate the $`Q^2`$ contribution. The constraint increases in complexity as additional terms are included. A comparison between the Kubon parametrization and our LGE fit with $`\lambda _M=2`$ is shown in Fig. 20. The Kubon analysis included only the data indicated by filled circles and was limited to $`Q^2<4`$ (GeV/$`c`$)2 while our analysis used all data shown and extended to 10 (GeV/$`c`$)2. The LGE parametrization fits the data well over a broader range and is compatible with pQCD, whereas the continued-fraction parametrization behaves badly soon after the range fitted by Kubon et al. By not imposing the pQCD constraint, the extrapolation deteriorates quite quickly. We believe that the continued-fraction method also underestimates the uncertainty in the rms radius due to the strong correlations among its parameters and its built-in normalization constraint at $`Q^2=0`$.

A rather different phenomenological parametrization can be made in the context of the vector meson dominance model at modest $`Q^2`$ matched to pQCD at large $`Q^2`$, denoted VMD+pQCD. That approach was pioneered by Gari and Krümpelmann  and recently refined by Lomon . Similarly, the classic dispersion-theory analysis of Höhler et al has recently been updated by Mergell et al to handle better the requirements of unitarity and the approach to the pQCD limit. These approaches have the advantage that all four electromagnetic form factors are analyzed simultaneously, thereby relating their isospin structure to an underlying model, and can be extended to timelike momentum transfer . By contrast, our approach is limited to spacelike momentum transfer and must construct the isospin form factors from four independent fits to individual form factors. Both the VMD and dispersion theory approaches appear to be capable of fitting the data as well as our linear expansion analysis, although the data have improved considerably since the analysis of Ref. . However, we omit detailed comparisons here because these models do not consider radial densities.

### F Importance of $`G_{En}`$ data at higher $`Q^2`$

The present data for $`G_{En}`$ do not extend high enough in $`Q^2`$ to determine the interior charge density as accurately for the neutron as for the proton or to permit reliable extrapolation to the scaling regime, but new data expected from an approved proposal  at Jefferson Laboratory should help considerably. The impact of extending the $`Q^2`$ range to 3.4 (GeV/$`c`$)2 is illustrated in Fig. 21. This analysis was performed using $`\lambda _E=0`$, which permits the greatest latitude at high $`Q^2`$. The left column shows the form factor and density fitted to published data (shown by open symbols), while the middle and right columns show the effect of pseudodata (shown as filled symbols) for two hypothetical scenarios. The middle scenario assumes that the new data would follow the Galster parametrization while the right scenario assumes that the new data would fall more rapidly than the dipole form factor for $`Q^2>2`$ (GeV/$`c`$)2. Both scenarios are compatible with the uncertainties extrapolated from the fit to the present data and can be fitted well, with error bands for $`Q^24`$ (GeV/$`c`$)2 consistent with the anticipated experimental precision. The reduction of the incompleteness error obtained by extending the measurements to higher $`Q^2`$ greatly improves the precision of the interior charge density. If measurements at higher $`Q^2`$ come close to the Galster parametrization, then the error band would be reduced in width with little change in its centroid. On the other hand, if new measurements of $`G_{En}/G_D`$ decrease with $`Q^2`$ in a manner similar to the proton charge form factor, the softer charge density would be reduced in the interior, moving toward the lower edge of the present error band. These scenarios have quite different asymptotic values for $`Q^4G_{En}`$, but the present data cannot distinguish between them. Furthermore, the recent VMD+pQCD analysis by Lomon  suggests that $`G_{En}/G_D`$ could reach an asymptotic value substantially higher than predicted by the Galster parametrization. Therefore, it is very important to extend $`G_{En}`$ data as far as possible in $`Q^2`$.

## VI Summary and Conclusions

We have employed expansions in complete sets of radial basis functions to parametrize nucleon Sachs form factors in terms of charge and magnetization densities. Our selection of data emphasizes recent polarization data. The inversion from form factor to density is based upon relativistic models in which the spatial frequency $`k=Q/\sqrt{1+\tau }`$ in the rest frame is related to the momentum transfer $`Q`$ in the Breit frame by Lorentz contraction. The maximum possible frequency sampled by electron scattering is then $`k<k_m`$ where $`k_m=2m`$ is determined by the nucleon Compton wavelength. A variety of models produce inversion formulas of the form

$$\stackrel{~}{\rho }(k)=G(Q^2)(1+\tau )^\lambda $$

but differ in the choice of $`\lambda `$. By considering the asymptotic form of $`G(Q^2)`$, we can limit the exponent $`\lambda `$ to 0,1,2. The relativistic soliton model suggests $`\{\lambda _E=0,\lambda _M=1\}`$, the original quark cluster model suggested $`\lambda _E=\lambda _M=1`$, and a more symmetric version of the quark cluster model gives choices $`\lambda _E=\lambda _M=2`$ that are compatible with pQCD without the somewhat artificial constraints upon $`\stackrel{~}{\rho }(k_m)`$ needed by the other models. In most of this paper we parametrized the Sachs form factors using the Laguerre-Gaussian expansion (LGE) and derived densities using $`\lambda _E=\lambda _M=2`$, but we have also analyzed the impact of the discrete ambiguity in $`\lambda `$ upon the radial densities. Although some of the details of the radial densities are affected by the discrete ambiguity in the relativistic inversion formula, their qualitative features are independent of $`\lambda `$.

We find that virtually identical fits to the Sachs form factors are obtained with either LGE or Fourier-Bessel expansions (FBE) and within a wide range these fits are independent of details of the parametrization, such as number of terms, radial scale parameters, and tail bias. The fitting procedure uses $`\stackrel{~}{\rho }(k)`$ pseudodata for $`k>k_{max}`$ to estimate the incompleteness error in radial density due to the limitation of experimental data for the range $`k<k_{max}<k_m`$. For a given choice of $`\lambda `$, the radial densities fitted with either LGE or FBE expansions are practically identical and are insensitive to details of the analysis. We find that the proton charge density is significantly broader than its magnetization density, consistent with the observation in recent recoil-polarization measurements that $`G_{Ep}/G_{Mp}`$ decreases in an almost linear fashion for $`1<Q^2<6`$ (GeV/$`c`$)2. Our result for the proton Sachs charge radius is consistent with a recent determination based upon the 1S Lamb shift in hydrogen. Similarly, we find the magnetization density is slightly broader for the neutron than for the proton. Each of these three densities exhibits a secondary peak in $`r^2\rho (r)`$ near $`11.5`$ fm that cannot be suppressed without sacrificing the fit to the corresponding form factor for $`Q^2>1`$ (GeV/$`c`$)2. This structure appears somewhat stronger for $`\lambda =2`$ than for $`\lambda =0`$, but $`\lambda =2`$ provides the clearest extrapolation to the pQCD limit.

The recent recoil-polarization data for $`F_{2p}/F_{1p}`$ appear to favor scaling with $`Q^1`$ rather than the $`Q^2`$ expected from quark helicity conservation in pQCD. Although that observation has stimulated some speculation about violations of quark helicity conservation due to orbital angular momentum or imposition of Poincaré invariance, in our analysis with $`\lambda _E=\lambda _M=2`$ we find that although $`QF_{2p}/F_{1p}`$ appears to be nearly constant for $`2<Q^2<6`$ (GeV/$`c`$)2 we nevertheless obtain a constant asymptotic value for $`Q^2F_{2p}/F_{1p}`$ for $`Q^220`$ (GeV/$`c`$)2 where $`G_{Mp}`$ scales with $`Q^4`$. Therefore, we find that the data are consistent with a broad maximum in $`QF_{2p}/F_{1p}`$ and do not require true $`Q^1`$ scaling.

We have compared the LGE parametrization for $`G_{En}`$ to fits based upon the Galster parametrization. Although the range of $`Q^2`$ remains too small to discriminate between these models, the Galster parametrization cannot be inverted using a relativistic relationship between intrinsic spatial frequency and Breit-frame momentum transfer unless $`\lambda 1`$. The traditional nonrelativistic inversion of the Galster form factor produces a charge density with an unphysical cusp at the origin while the relativistic fit to the data using the LGE form factor softens the interior density and removes the cusp. However, the incompleteness error in the neutron charge density remains fairly large because the available $`G_{En}`$ data are limited to small $`Q^2`$ and the data for $`0.4<Q^2<1.6`$ (GeV/$`c`$)2 have relatively large uncertainties. More precise data for $`0.5<Q^2<1.5`$ (GeV/$`c`$)2 are expected soon and an experiment for $`Q^23.4`$ (GeV/$`c`$)2 is in preparation. These data should improve the accuracy of the neutron charge density considerably, but data approaching 20 (GeV/$`c`$)2 will probably be needed to test scaling in the neutron and in the isospin form factors.

Combining the neutron and proton charge densities, we deduced the up and down quark radial distributions assuming isospin symmetry and neglecting heavier quarks. This schematic model suggests that the distribution is slightly broader for up quarks than for down quarks in the proton. With $`\lambda _E=2`$ we also observe a statistically significant negative density for down quarks near 1 fm that might be attributed to the $`\overline{d}`$ content of the pion cloud.

Although we cannot claim that there is a unique relationship between form factors and densities, expansion of densities in a complete radial bases provides physically appealing parametrizations of form factor data that are applicable over a wide range of $`Q^2`$. The use of linear expansions in complete bases minimizes the model dependence of the fitted form factors and provides more realistic error bands in both spatial and momentum representations. Therefore, even if the identification with static densities is discounted, the fitted densities do provide useful parametrizations of the form factors nonetheless. The choice $`\lambda _E=\lambda _M=2`$ automatically satisfies pQCD scaling and provides a natural means for extrapolating form factors to higher $`Q^2`$ for the purpose of planning future experiments.

###### Acknowledgements.

We thank O. Gayou for providing a table of $`G_{Ep}`$ data and X. Ji, C. Perdrisat, R. Madey, J. Friar, and R. Rosenfelder for useful discussions. The support of the U.S. National Science Foundation under grant PHY-9971819 is gratefully acknowledged.