Measurement of the neutron charge radius and the role of its constituents

The neutron is a cornerstone in our depiction of the visible universe. Despite the neutron zero-net electric charge, the asymmetric distribution of the positively- (up) and negatively-charged (down) quarks, a result of the complex quark-gluon dynamics, lead to a negative value for its squared charge radius, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\langle {r}_{{\rm{n}}}^{2}\rangle$$\end{document}⟨rn2⟩. The precise measurement of the neutron’s charge radius thus emerges as an essential part of unraveling its structure. Here we report on a \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\langle {r}_{{\rm{n}}}^{2}\rangle$$\end{document}⟨rn2⟩ measurement, based on the extraction of the neutron electric form factor, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${G}_{{\rm{E}}}^{{\rm{n}}}$$\end{document}GEn, at low four-momentum transfer squared (Q2) by exploiting the long known connection between the N → Δ quadrupole transitions and the neutron electric form factor. Our result, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\langle {r}_{{\rm{n}}}^{2}\rangle =-0.110\pm 0.008\,({{\rm{fm}}}^{2})$$\end{document}⟨rn2⟩=−0.110±0.008(fm2), addresses long standing unresolved discrepancies in the \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\langle {r}_{{\rm{n}}}^{2}\rangle$$\end{document}⟨rn2⟩ determination. The dynamics of the strong nuclear force can be viewed through the precise picture of the neutron’s constituent distributions that result into the non-zero \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\langle {r}_{{\rm{n}}}^{2}\rangle$$\end{document}⟨rn2⟩ value.

T he study of the nucleon charge radius has been historically instrumental towards the understanding of the nucleon structure. In the neutron case, it is the highly complicated dynamics of the strong force between quarks and gluons, the fermionic nature of quarks and spin-orbit correlations that leads to an asymmetric distribution of u-and d-quarks in it, thus resulting in a negative value for hr 2 n i. The precise measurement of hr 2 n i becomes a critical part of our understanding of the nucleon dynamics. Furthermore, employing new, different techniques in extracting this fundamental quantity has proven most valuable, as recently exhibited in the proton's case: the disagreement of the proton charge radius, r p , as determined using the Lamb shift measurement in the muonic hydrogen atom 1 , with the earlier results based on the hydrogen atom and the electron scattering measurement, gave rise to the proton radius puzzle 2 . In turn, this led to a significant reassessment of the methods and analyses utilized in the proton radius extraction, and to the consideration of physics beyond the standard model as potential solutions to this discrepancy. Various atomic and nuclear physics techniques were employed for the proton r p measurement. However, in the neutron case, the determination of hr 2 n i is more challenging since no atomic method is possible and the electron scattering method suffers from severe limitations due to the absence of a free neutron target. Thus, the extraction of hr 2 n i has been uniquely based on the measurement of the neutron-electron scattering length b ne , where low-energy neutrons are scattered by electrons bound in diamagnetic atoms. The hr 2 n i measurements adopted by the particle data group (PDG) 3-6 exhibit discrepancies, with the values ranging from hr 2 n i ¼ À0:114 ± 0:003 4 to hr 2 n i ¼ À0:134 ± 0:009 ðfm 2 Þ 5 . Among the plausible explanations that have been suggested for this, one can find the effect of resonance corrections and of the electric polarizability, as discussed e.g., in ref. 4 . However, these discrepancies have not been fully resolved, a direct indication of the limitations of this method.
An alternative way to determine hr 2 n i is offered by measuring the slope of the neutron electric form factor, G n E , at Q 2 → 0, which is proportional to hr 2 n i. In the past, determinations of G n E at finite Q 2 were typically carried out by measuring double polarization observables in quasi-elastic electron scattering from polarized deuterium or 3 He targets using polarized electron beams [7][8][9][10][11][12][13][14][15][16][17][18][19][20][21] . However, these measurements were not able to access G n E at a sufficiently low Q 2 range so that the slope, and subsequently the hr 2 n i can be determined. In this work we rely on an alternative path to access G n E . It has long been known 22,23 that the ratios of the quadrupole to the magnetic dipole transition form factors (TFFs) of the proton, C2/M1 (CMR) and E2/M1 (EMR), are related to the neutron elastic form factors ratio G n E =G n M . Here, we follow that path and we access G n E at low momentum transfers from high precision measurements of the two quadrupole TFFs. The main steps of this work are summarized here for clarity. First, we extract G n E from the quadrupole TFF data, at low momentum transfers [24][25][26][27][28] , utilizing the form factor relations 22,23 determined within the SU(6) and the large-N c frameworks. The variance of the G n E results from the two analyses is treated as a theoretical uncertainty. The G n E ðQ 2 Þ form factor is then parametrized and fitted to the data, and hr 2 n i is determined from the G n E -slope at Q 2 = 0. Finally, we perform the flavor decomposition of the neutron and the proton form factors measurements and derive the flavor dependent quark densities in the nucleon, which reveal with high precision the role of the quark contributions to the neutron charge radius.

Results
A consequence of the SU(6) spin and flavor symmetry group that relates the nucleon and the Δ resonance leads to the following expression 22 where |q| is the virtual photon three-momentum transfer magnitude in the γN center of mass frame and M N is the nucleon mass. The n b parametrizes the contribution from threequark current terms, that tend to slightly increase the C2/M1 ratio (or correspondingly decrease the G n E =G n M ), an SU(6) symmetry breaking correction that has been theoretically quantified to~10% 22 (i.e., n b~1 .1). If one chooses to follow the most conservative path, a theoretical uncertainty can be assigned to this term that is equal to the full magnitude of the symmetry breaking contributions i.e., n b = 1.1 ± 0.1. Considering the confidence with which the underlying theory is able to determine the level of the symmetry breaking contributions, the above assumption leads to a safe estimation, and most likely to an overestimation, of the theoretical uncertainty.
In one of the first steps of this work we check the validity of the underlying theory using experimental measurements. The wealth of the TFF 24-32 and of the G n E =G n M 7-21 world data allow to quantify the magnitude of the symmetry breaking corrections from the analysis of the experimental measurements. In Fig. 1a we show the neutron G n E =G n M world data [7][8][9][10][11][12][13][14][16][17][18][19][20][21] (open black cirles), and we compare it to the G n E =G n M ratios that we have  33 are the same as in panel (a). The error bars correspond to the total uncertainty, at the 1σ or 68% confidence level. ARTICLE NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-021-22028-z derived from the TFFs C2/M1 measurements 24-32 (filled boxes) utilizing Eq. (1) with n b = 1 (i.e., uncorrected for the symmetry breaking contributions). By parametrizing the two data sets and then forming their ratio we can experimentally determine the magnitude of the n b (Q 2 ) contribution. A variety of functional forms have been explored to identify the functions that can provide a good fit to the data. All the appropriate functions that offer a good fit have been considered in the determination of n b and the variance of the results arising from the choice of the functional form is adopted as an uncertainty. The procedure is further refined using lattice Quantum Chromodynamics (LQCD) results at low momentum transfers, where neutron data do not exist. In particular, we extracted the ratio G n E =G n M from numerical simulations within LQCD using the G n E and G n M data of ref. 33 . The LQCD data provide further guidance on the Q 2 -dependence of the G n E =G n M ratio based on ab-initio QCD calculations, in a region where neutron form factor data are not available. The LQCD input results to a rather small refinement of ≤0.003 in the determination of n b . The details on the determination of n b (Q 2 ) are given in Section 2.1 of the Supplementary Information. The experimentally derived n b (Q 2 ) is found in excellent agreement with the theoretical prediction 22 as seen in Fig. 2; this in-turn offers further credence to the theoretical effort in ref. 22 . Furthermore, the fitted parametrizations allow to constrain the n b (Q 2 ) uncertainty by a factor of two compared to the most conservative n b = 1.1 ± 0.1 (i.e., as indicated by the width of the uncertainty band in Fig. 2), but also to determine these contributions accurately at very low momentum transfers where the analysis of the current TFF data takes place for the G n E extraction. The LQCD results entering our analysis are compared to the experimental world data and they exhibit a very good agreement as shown in Fig. 1a. The parameters of the LQCD calculation are such that they reproduce the physical value of the pion mass. Thus, such a calculation eliminates a major source of systematic uncertainties, that is, the need of a chiral extrapolation. Furthermore, the lattice results include both the connected and disconnected diagrams, and therefore G n E and G n M include both valence and sea quark contributions.
In our analysis we have extracted the G n E under two scenarios; in one case we consider the conservative path where n b = 1.1 ± 0.1, while in the second we utilize the n b as we have determined it from the experimental world data. The two sets of results come to an agreement at the ≤3% level; this is much smaller than to the overall G n E uncertainty. A slightly improved G n E uncertainty is obtained in the latter case due to an improved level of the n b uncertainty, when these contributions are determined from the world data. The G n E uncertainty of our results is driven by the following sources: (i) Experimental (statistical and systematic) uncertainties in the determination of the C2/M1 ratio. (ii) Uncertainties in the determination of C2/M1 due to the presence of non-resonant pion electro-production amplitudes that interfere with the extraction of the resonant amplitudes. These effects were studied by employing theoretical pion electro-production models in the data analysis; they were further investigated experimentally by measuring C2/M1 through an alternative reaction channel, the pðe; e 0 pÞγ 28 , where one employs a different theoretical framework for the ratio extraction (see Supplementary Information, Section 1). (iii) The uncertainty of the symmetry breaking terms δn b , as discussed above. (iv) The uncertainty introduced by the choice of the G n M -parametrization, in order to extract the G n E from the G n E =G n M ratio (as typically done in such cases e.g., 9,17 ). In this work we have used the one from ref. 34 and we have quantified the associated uncertainty by repeating the analysis with alternative parametrizations. We have found a~0.5% effect, which is rather small compared with the total G n E uncertainty. The G n E results are displayed in Fig. 3a.  The relation between G n E and the quadrupole transition form factors has also been established through large-N c relations 23 . The relations take the form where F ðpnÞ 2 are the nucleon Pauli form factors, M Δ is the mass of the Δ, and Here one is free from any additional correction terms, such as the symmetry breaking contributions of Eq. (1). Another advantage is that the experimental database is extended to include the Electric quadrupole (E2) transition, which in turn allows for an improved extraction of G n E . Being able to extract G n E independently through the Coulomb and the Electric quadrupole transitions offers a strong experimental test to the validity of the large-N c relations and allows to quantify their level of theoretical uncertainty. The above relations come with a 15% theoretical uncertainty 23 that is treated accordingly in the G n E analysis. The G n E extraction from the Coulomb and from the Electric quadrupole transitions agree nicely within that level, as can be seen in Fig. 1b 34 . For G p E we performed an updated parametrization so that we may include recent measurements from ref. 35 that were not yet available in ref. 34 (see Supplementary Information, Section 4). For the large-N c analysis the final results integrate both of the quadrupole transition form factors from each experiment, when both of them were simultaneously measured, into one G n E measurement (see Supplementary  Information, Section 2.2). The extracted G n E results from the large-N c analysis are displayed in Fig. 3a and are compared to the results from the SU(6) analysis in the same figure. The SU(6) analysis is in agreement with the large-N c analysis G n E results. For our final G n E result we consider the weighted average of the two values, as shown in Fig. 3b. The variance of the two values is treated as an additional G n E theoretical uncertainty, and is accounted for accordingly in the hr 2 n i extraction. The neutron mean square charge radius is related to the slope of the neutron electric form factor as Q 2 → 0 through In order to determine the charge radius the data have to be fitted to a functional form, and the slope has to be determined at Q 2 = 0. It is important that a proper functional form is identified so that model dependent biases to the fit are avoided. In the past, the experimental data would allow to explore functional forms for G n E ðQ 2 Þ with only two free parameters, that were lacking the ability to determine the neutron charge radius. The updated data allow to introduce an additional free parameter and to extract the hr 2 n i from measurements of the G n E that was not possible previously. Our studies have shown that is the most robust function for the radius extraction, where τ ¼ Q 2 =4M 2 N , and A, B, C are free parameters (see Supplementary  Information, Section 3). Our fits employ the G n E data discussed in this work as well as the G n E world data from [7][8][9][10][11][12][13][14][15][16][17][18][19][20][21] . The function describes the data very well, with a reduced χ 2 of 0.74. The parameters obtained are A = 0.505 ± 0.079 (GeV/c) 2 , B = 1.655 ± 0.126, C = 0.909 ± 0.583, and Q 2 in units of (GeV/c) 2 , leading to a value of hr 2 n i ¼ À0:110 ± 0:008 ðfm 2 Þ, as shown in Fig. 4. When the uncertainty of the symmetry breaking contributions in the SU(6) analysis is treated conservatively (i.e., n b = 1.1 ± 0.1) the final result becomes hr 2 n i ¼ À0:109 ± 0:009 ðfm 2 Þ with a reduced χ 2 of 0.74. Here we observe that the hr 2 n i-uncertainty is not affected significantly by the different treatment of the symmetry breaking contributions in the two cases.
The charge radius extraction is further explored through fits that are constrained within a limited range at low Q 2 where G n E remains monotonic, namely from Q 2 = 0 to 0.4 (GeV/c) 2 . In the fits to the data, the functional forms can be divided into two groups, those based on polynomials with varying orders and those that are based on rational forms (see Supplementary Information, Section 3.1). For the charge radius, the weighted average is extracted separately for each one of the two groups. A systematic uncertainty is also quantified within each group (i.e., a model uncertainty of the group) from the weighted variance of the results from all the fits within the group. The results from the two groups tend to have a similar overall uncertainty. A small systematic difference of the two group's central hr 2 n i values is observed, as studies over a varying fitting range have shown. For that reason a third uncertainty is determined: here we consider the spread of the two central values as indicative of the uncertainty that is associated with the choice of the group. Therefore, the final result is given by the average of the two group values for hr 2 n i, while the half of the difference of the two values is assigned as an additional uncertainty. The details of the studies are presented in the Supplementary Information, Sections 3.1 and 3

.2.
The results from the low-Q 2 fits for all groups of functions and for variations to the fitting range are shown in the Supplementary Information, Table 7 and Table 8. The low-Q 2 fits confirm the hr 2 n i extraction from the fits over the complete G n E database, but they are vulnerable to model uncertainties that are associated with the choice of the fitted parametrization and they are not able to improve the hr 2 n i extraction. The studies indicate that more G n E measurements are needed at lower momentum transfers so that a more competitive extraction can become possible from the low-Q 2 fits. This comes as no surprise when one considers the corresponding case for the proton, in which case the charge radius extraction from fits within a limited Q 2 -range required measurements at significantly smaller momentum transfers, namely at Q 2 = 0.0002 (GeV/c) 2 − 0.06 (GeV/c) 2 35 . Fig. 4 The neutron mean square charge radius. The hr 2 n i measurement from this work (red circle), with the error bar corresponding to the total uncertainty at the 1σ or 68% confidence level, and from references 3-6 (black box) included in the PDG analysis for hr 2 n i. The orange-band indicates the PDG averaged hr 2 n i value. The new weighted average of the world data is also shown (blue diamond) when the new hr 2 n i measurement reported in this work is included in the calculation. ARTICLE NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-021-22028-z The quadrupole TFF data at low momentum transfers [24][25][26][27][28] , acquired at MAMI and at JLab/Hall-A, provide the critical G n E data that were missing lower than Q 2 = 0.20 (GeV/c) 2 and make the hr 2 n i extraction possible. We now explore the potential of extending the current analysis to higher momentum transfers. This study aims to observe the effect of the hr 2 n i extraction if we enrich the G n E database with measurements in the region where G n E data already exist, higher than Q 2 = 0.20 (GeV/c) 2 . The relations between the G n E and the quadrupole transition form factors hold on very solid ground in the low Q 2 region, i.e., lower than Q 2 = 0.20 (GeV/c) 2 . On the other hand, they tend to hold less well at high momentum transfers. The relations do not come with a sharp Q 2 cut-off-value after which they do not hold, but one should avoid the Q 2 = 1 (GeV/c) 2 region where larger theoretical uncertainties could bias the charge radius extraction. In extending the database higher in Q 2 , care has to be given so that any additional data at intermediate momentum transfers will benefit the fits without compromising the hr 2 n i extraction by the gradually increasing theoretical uncertainties. Our studies showed that when we integrate in the analysis the G n E data that we extract from the CLAS measurements up to Q 2 = 0.52 (GeV/c) 2 31 (see Fig. 3b, green boxes) we find that hr 2 n i ¼ À0:107 ± 0:007 ðfm 2 Þ, compared with hr 2 n i ¼ À0:110 ± 0:008 ðfm 2 Þ when these additional data are not included (see the Supplementary Information, Section 3.2 for details). Last, when the same data-set is included in the fits within the limited low-Q 2 range, as discussed in the previous paragraph, we find that hr 2 n i ¼ À0:111 ± 0:006 ± 0:002 mod ± 0:004 group ðfm 2 Þ. Here the last two uncertainties (mod and group) are model-related uncertainties associated with the choice of the fitted parametrization (see Supplementary Information, Section 3.2 for details). We do not observe any additional benefit by extending the measurements higher in Q 2 since the fits' uncertainties do not improve when more data are included up to Q 2 = 1 (GeV/c) 2 . In conclusion, a small benefit to the charge radius uncertainty can be observed when additional data up to Q 2 = 0.5 (GeV/c) 2 are utilized for the charge radius extraction. If one decides to eliminate any risk of introducing theoretical bias from the inclusion of the intermediate momentum transfer measurements for the final result, one can conservatively adopt the analysis that does not include these additional data, namely hr 2 n i ¼ À0:110 ± 0:008 ðfm 2 Þ.

Discussion
Our analysis and results offers valuable input towards addressing long standing unresolved hr 2 n i discrepancies of the b ne -measurements, which display a ≈10% tension between the results, suggesting that there are still unidentified systematic uncertainties associated with this method of extraction. Our measurement is in disagreement with ref. 5 and supports the results of refs. 3,4 . Considering that here we cross check hr 2 n i using a different extraction method, there is a strong argument so as to exclude the value of ref. 5 from the world data average. In such a case, the new weighted average value of the world data when we include our measurement and we exclude the one of ref. 5 , becomes hr 2 n i ¼ À0:1152 ± 0:0017 ðfm 2 Þ. Based on the current work, the particle data book value of hr 2 n i ¼ À0:1161 ± 0:0022 ðfm 2 Þ is adjusted by~1% and improves its uncertainty by~23%. We also note that our result agrees very well with a recent hr 2 n i calculation that is based on the determination of the deuteron structure radius in chiral effective field theory and utilizes atomic data for the difference of the deuteron and proton charge radii 36 .
The neutron's non-zero mean charge radius is a direct consequence of the asymmetric distribution of the positively-charged (up) and of the negatively-charged (down) quarks in the system, a consequence of the non-trivial quark gluon dynamics of the strong force. The quark distributions offer a detailed view as to how the non-zero hr 2 n i value arises. Here, one has to work on the infinite-momentum frame 37 since it offers the inherent advantage that a true transverse charge density can be properly defined as the matrix element of a density operator between identical initial and final states. We find that the results of our analysis on hr 2 n i are particularly sensitive to the neutron's long-distance structure, and offer a significant improvement (factor of 2) in the precision of the neutron charge density at its surface (see Supplementary  Information Fig. 9). We extract the neutron and the proton charge densities at the infinite-momentum-frame from the most recent nucleon form factor parametrizations, where for G n E we use the one determined in this work. The details are presented in the Supplementary Information, Section 4. The extracted neutron and proton charge densities are shown in Fig. 5. From the two nucleon densities, invoking charge symmetry, and neglecting the s s contribution, we derive the uand d-quark densities with an improved precision as shown in Fig. 5 (see Supplementary Information, Section 4 for details). The flavor dependent densities show that the singly-represented quark in the nucleon has a wider distribution compared with the doubly-represented quarks, which in turn exhibit a larger central quark density. Although the concentration of the two negatively charged quarks at the center of the neutron may appear to contradict the negative sign of the neutron's mean square charge radius, this is not truly the case. The 3D Breit frame and 2D infinite-momentum distributions are directly related to each other and the apparent discrepancies between the distributions in the two frames simply result from kinematical artifacts associated with spin 38 . The effect is rather dramatic in the neutron, where the rest-frame magnetization is large and negative. The contribution it induces competes with the convection contribution and gradually changes the sign at the center of the charge distribution as one increases the momentum of the neutron. Thus, the appearance of a negative region around the center of the neutron charge distribution in the infinitemomentum frame is just a manifestation of the contribution induced by the rest-frame magnetization.
In conclusion, we report on an alternative measurement of the neutron charge radius, based on the measurement of the neutron electric form factor G n E . An alternative path to the measurements based on the scattering of neutrons by electrons bound in diamagnetic atoms is presented. Our value of hr 2 n i ¼ À0:110 ± 0:008 ðfm 2 Þ offers valuable input towards addressing long standing unresolved discrepancies in the hr 2 n i measurements, rejects earlier measurements, and improves the precision of the hr 2 n i world data average value. Furthermore, our data offer access to the associated dynamics of the strong nuclear force through the precise mapping of the quark distributions in the neutron that contribute to its non-zero charge radius. The current work lays the path for hr 2 n i measurements of higher precision. New experimental proposals based on this method, e.g., Jefferson Lab LOI 12-20-002, offer to improve the precision of the hr 2 n i measurement by nearly a factor of 2. Future experimental efforts will be able to utilize upgraded experimental setups that will fully exploit the advantages of this method. In particular, pushing the low momentum transfer limits of high precision measurements can lead to a further improvement in the precision of the hr 2 n i extraction.

Methods
Here we extract G n E , at low momentum transfers, from measurements of the quadrupole transition form factors. The neutron charge radius is then extracted from the slope of G n E at Q 2 = 0. Utilizing the G n E data and the world data for the nucleon elastic form factors the flavor decomposition of the nucleon electromagnetic form factors is performed, and the u-and d-quark distributions in the nucleon are extracted. The main steps of this work are as follows: 1. We extract G n E from the Coulomb quadrupole and the Electric quadrupole transition form factor data at low momentum transfers [24][25][26][27] utilizing the form factor relations (Eqs. (1-3)) that have been determined within the SU (6) 22  3. The neutron mean square charge radius hr 2 n i is obtained from Eq. (4) by fitting the G n E data to the functional form in Eq. (5) and determining the slope at Q 2 = 0. This functional form was shown to be the most robust function for the radius extraction from the neutron data. The fit employs additional G n E data reported here and the G n E world data extending to higher momentum transfers. 4. The neutron and proton densities are derived at the infinite momentum frame through ρðbÞ ¼ where b is the transverse distance, τ = Q 2 /4m 2 and J 0 the 0 th order cylindrical Bessel function. Here we utilize the most recent parametrizations for the nucleon form factors, where for G n E we use the one derived in this work. From the neutron and proton densities, invoking charge symmetry, and neglecting the s s contribution, we then extract the uand d-quark densities in the proton (or doubly-represented and singly-represented quarks in the nucleon, respectively), where ρ u ðbÞ ¼ ρ p ðbÞ þ ρ n ðbÞ=2 ð7Þ and ρ d ðbÞ ¼ ρ p ðbÞ þ 2ρ n ðbÞ: ð8Þ

Data availability
All the relevant data in this work are available from the authors upon request. The data for the quadrupole TFFs used in this work are publicly available in their original publications [24][25][26][27][28]31 where they are described in detail. The raw data from these experiments are archived in Jefferson Laboratory's mass storage silo and at Temple University, Department of Physics. The G n E data that we have used for the hr 2 n i extraction are available in the following sources: (i) The new G n E data that we have derived in this work, from the analysis of the quadrupole TFFs, are available in the Supplementary  Information document. (ii) The previously published G n E world-data are available in their original publications [7][8][9][10][11][12][13][14][15][16][17][18][19][20][21] .