Electromagnetic charge radius of the pion at high precision

TL;DR

This paper determines the pion charge radius with high precision using a novel formalism based on analyticity and unitarity, combining timelike and spacelike data without relying on specific parametrizations.

Contribution

It introduces a new method that leverages phase information and high-precision data to accurately determine the pion charge radius, improving previous precision significantly.

Findings

Charge radius of the pion is $0.657 m$ with an uncertainty of $0.003 m$

Precision improved by a factor of about 2.7 over previous averages

Method does not depend on the unknown phase above the inelastic threshold

Abstract

We present a determination of the pion charge radius from high precision data on the pion vector form factor from both timelike and spacelike regions, using a novel formalism based on analyticity and unitarity. At low energies, instead of the poorly known modulus of the form factor, we use its phase, known with high accuracy from Roy equations for $\pi\pi$ elastic scattering via the Fermi-Watson theorem. We use also the values of the modulus at several higher timelike energies, where the data from $e^+e^-$-annihilation and $\tau$-decay are mutually consistent, as well as the most recent measurements at spacelike momenta. The experimental uncertainties are implemented by Monte-Carlo simulations. The results, which do not rely on a specific parametrization, are optimal for the given input information and do not depend on the unknown phase of the form factor above the first inelastic threshold. Our prediction for the charge radius of the pion is $r_\pi=(0.657 \pm 0.003) \fm $, which amounts to an increase in precision by a factor of about 2.7 compared to the PDG average.

Tables

Table 1. Table 1: Central values and standard deviations for the quantity ⟨ r π 2 ⟩ delimited-⟨⟩ superscript subscript 𝑟 𝜋 2 \langle r_{\pi}^{2}\rangle obtained by the averaging prescription described in the text for each experiment.

	Bern phase	Madrid phase
CMD2 06	$0.4281\pm 0.0064$	$0.4279\pm 0.0061$
SND 06	$0.4323\pm 0.0059$	$0.4327\pm 0.0055$
BABAR 09	$0.4343\pm 0.0048$	$0.4351\pm 0.0046$
KLOE 11	$0.4304\pm 0.0055$	$0.4304\pm 0.0048$
KLOE 13	$0.4311\pm 0.0060$	$0.4313\pm 0.0054$
BESIII 15	$0.4293\pm 0.0063$	$0.4319\pm 0.0057$
CLEO 00	$0.4340\pm 0.0060$	$0.4346\pm 0.0054$
ALEPH 05	$0.4315\pm 0.0067$	$0.4318\pm 0.0064$
Belle 05	$0.4347\pm 0.0056$	$0.4356\pm 0.0051$
OPAL 12	$0.4266\pm 0.0082$	$0.4265\pm 0.0079$

Equations

$$r_{\pi }=(0.672\pm 0.008)\text{fm},$$

$$\frac{1}{\pi }{\int }_{t_{\text{in}}}^{\infty }\frac{dt}{t}|F(t){|}^2\le I,$$

$$\tilde{z}(t)=\frac{\sqrt{t_{\text{in}}}-\sqrt{t_{\text{in}}-t}}{\sqrt{t_{\text{in}}}+\sqrt{t_{\text{in}}-t}},$$

$$g(z)=F(\tilde{t}(z))[\mathcal{O}(\tilde{t}(z)){]}^{-1}w(z)\omega (z).$$

$$\mathcal{O}(t)=\exp \left(\frac{t}{\pi }{\int }_{4m_{\pi }^2}^{\infty }d{t}^{\prime }\frac{\delta (t^{\prime })}{t^{\prime }(t^{\prime }-t)}\right),$$

$$w(z)=\sqrt{\frac{1-z}{1+z}},$$

$$\omega (z)=\exp \left(\frac{\sqrt{t_{\text{in}}-\tilde{t}(z)}}{\pi }\underset{t_{\text{in}}}{\overset{\infty }{\int }}\frac{\ln |\mathcal{O}(t^{\prime })|\mathrm{d}{t}^{\prime }}{\sqrt{t^{\prime }-t_{\text{in}}}(t^{\prime }-\tilde{t}(z))}\right).$$

$$\det\left|\begin{array}[]{c c c}\vspace{0.2cm}I-g(0)^{2}-g^{\prime}(0)^{2}&\bar{g}(z_{a})&\bar{g}(z_{b})\\ \vspace{0.2cm}\bar{g}(z_{a})&\displaystyle\frac{z^{4}_{a}}{1-z^{2}_{a}}&\displaystyle\frac{(z_{a}z_{b})^{2}}{1-z_{a}z_{b}}\\ \bar{g}(z_{b})&\displaystyle\frac{(z_{a}z_{b})^{2}}{1-z_{a}z_{b}}&\displaystyle\frac{(z_{b})^{4}}{1-z_{b}^{2}}\\ \end{array}\right|\geq 0,$$

$$A⟨r_{\pi }^2{⟩}^2+2B⟨r_{\pi }^2⟩+C\ge 0,A\le 0.$$

$$\bar{\sigma }={\left[\sum _{i,j=1}^n(C(f{)}^{-1}{)}_{ij}\right]}^{-1/2},$$

$$r_{\pi }=(0.657\pm 0.003)\text{fm}.$$

Full text

Introduction.— The electromagnetic charge radius of the pion is a fundamental observable of the strong interactions, with a long history dating from over half a century. As with similar observables like the proton radius or the pion-nucleon $\sigma$ term, its accurate determination is crucial for precision tests of the Standard Model at low energy, especially of Chiral Perturbation Theory (ChPT) and lattice QCD.

In a relativistic theory, the mean charge radius squared is related to the slope at $t=0$ of the pion electromagnetic form factor $F(t)$, according to the expansion $F(t)=1+\frac{1}{6}\langle r_{\pi}^{2}\rangle t+O(t^{2})$. The most recent value of the pion charge radius, $r_{\pi}\equiv\langle r_{\pi}^{2}\rangle^{1/2}$, quoted by PDG [1]

[TABLE]

is obtained mainly from values of the form factor measured at small spacelike momenta $t<0$ from $e\pi\to e\pi$ and $eN\to e\pi N$ processes, extrapolated to $t=0$ using simple dipole-like parametrizations. ChPT at two-loop order [2] in the timelike region, $t>0$, yields the best value $\langle r_{\pi}^{2}\rangle=(0.437\pm 0.016)\,\text{fm}^{2}$, i.e. $r_{\pi}=(0.661\pm 0.015)\,\text{fm}$. The calculations in lattice QCD, summarized in [3], are consistent with these values, but have not yet reached the same precision.

There is a large amount of information on the pion form factor which is not used in the radius extractions quoted above. For $t>4m_{\pi}^{2}$, where $F(t)$ is a complex function, its modulus is measured either from $e^{+}e^{-}\to\pi^{+}\pi^{-}$ annihilation or, using isospin symmetry, from the $\tau^{-}\to\pi^{-}\pi^{0}\nu_{\tau}$ decay. This input involves however energies more distant from $t=0$, so its use for the radius extraction requires a nontrivial analytic continuation. The properties of the form factor that follow from causality and unitarity play here an important role. It is known that $F(t)$ is analytic in the complex $t$ plane with a unitarity cut along the region $(4m_{\pi}^{2},\infty)$ of the real axis, and its phase below the first inelastic threshold, in the limit of exact isospin symmetry, is equal by the Fermi-Watson theorem [4, 5] to the $P$-wave phase shift of $\pi\pi$ elastic scattering. Many analyses which exploit these properties have been performed over several decades, the pion form factor being actually one of the observables most investigated in dispersion theory. However, the standard dispersive relations, written only in terms of the phase (the so-called Omnès representation), or the modulus, or the imaginary part, always require some poorly known or unavailable input. Therefore, although in some papers (for instance Refs. [6, 7]) the reported uncertainty is quite small, a certain model dependence is unavoidable in the standard dispersive calculations of the pion radius.

In this work we apply an alternative mathematical formalism, proposed in [8] (see also [9]), which exploits in an optimal way analyticity, unitarity and the information available on the form factor. Essentially, the formalism is a mixed phase-modulus dispersive representation, which uses as input at low energies the phase, very precisely known from $\pi\pi$ scattering, and at higher energies the modulus, measured by high precision experiments. Using techniques of optimization theory for analytic functions, with no specific parametrization, we obtain for each input upper and lower bounds on the charge radius. The statistical distribution of the experimental input was then accounted for by Monte Carlo simulations, to convert the bounds into allowed intervals with definite confidence levels. This work considerably improves our previous study Ref. [10], by properly treating the experimental errors and by including more experimental data. Similar techniques have been applied recently in Ref. [11] for evaluating the low-energy hadronic contribution to muon $g-2$.

Theoretical and experimental inputs.— As shown in [12], the first significant inelastic threshold for the pion vector form factor is produced by the $\omega\pi$ channel at $\sqrt{t_{\text{in}}}=m_{\omega}+m_{\pi}=0.917\,\,\text{GeV}$. Below $t_{\text{in}}$ we take the phase of the form factor as the $P$-wave phase shift of $\pi\pi$ elastic scattering, which has been calculated very accurately from Roy equations and ChPT. We use the phase shift calculated in [13, 14] and the so-called “constrained” fit to data from [15], which we refer to as Bern and Madrid phases, respectively.

Above the inelastic threshold, where the Fermi-Watson theorem is no longer valid and the phase of the form factor is not known, we make use of the measurements of the modulus. In our conservative approach, we do not require the pointwise knowledge of $|F(t)|$, but adopt a weaker condition, expressed by a weighted integral of $|F(t)|^{2}$. Several weighted integrals have been investigated in the previous work [10], and stability of the results was proved for a large class of weights. Here we adopted the specific condition

[TABLE]

for which the available information allows an accurate and conservative estimate of the upper bound $I$. We evaluated the integral using the BABAR data [16] on the modulus $|F(t)|$ from $t_{\text{in}}$ up to $\sqrt{t}=3\,\,\text{GeV}$, smoothly continued with a constant value for the modulus in the range $3\,\,\text{GeV}\leq\sqrt{t}\leq 20\,\text{GeV}$, and a $1/t$ decreasing modulus at higher energies, as predicted by QCD scaling [17, 18]. This led to the value $I=0.578\pm 0.022$, where the uncertainty is due to the BABAR experimental errors. The contribution of the range above 3 GeV to the integral is actually of only $1\%$. As discussed in previous works [10, 11], the assumed behaviour of the modulus above 3 GeV largely overestimate perturbative QCD calculations [19], so the quoted value of $I$ is a conservative upper bound of the integral (2). This leads to weaker bounds on the charge radius, due to an exact monotony property of these bounds with respect to the value of $I$ [9, 10], which makes our predictions conservative. For the same reason, in the numerical analysis we have used as input for $I$ the quoted central value increased by the error.

The measurements of the modulus of $F(t)$ below the inelastic threshold are expected to further improve the precision. This input will be useful only if one can identify an energy range, sufficiently close to the point $t=0$ in order to play a role in the radius extraction, and where accurate measurements are available. Unfortunately, at very low energies the data have still large errors and there are discrepancies between the experiments. As shown in [10, 11], the above requirements are satisfied by the range $(0.65-0.71)\,\text{GeV}$. The number of experimental points in this range varies from experiment to experiment. For the $e^{+}e^{-}$-annihilation experiments, there are 2 points each from CMD2 [20] and SND [21], 26 points for BABAR [16, 22], 8 points each for KLOE 2011 [23] and KLOE 2013 [24] and 10 points for BESIII [25]. For the $\tau$-decay experiments, there are 3 points each for CLEO [26], ALEPH [28, 27] and OPAL [29], and 2 points for Belle [30]. We note that, compared to Ref. [10], we include now in addition the very recent data of BESIII [25], the alternative KLOE analysis [23], and the $\tau$-decay data from [26]-[30]. Several corrections, discussed in detail in Appendix B of Ref. [11], have been applied to the data in order to obtain the proper values of $|F(t)|$ in the isospin limit, required in the formalism.

As we mentioned, there is also rich experimental information on the pion form factor in the spacelike region, $t<0$. In our analysis we used only the values of the form factor measured recently with high-precision by Jefferson Lab $F_{\pi}$ Collaboration [31] at larger spacelike energies, $F(-1.60\,\,\text{GeV}^{2})=0.243\pm 0.012_{-0.008}^{+0.019}$ and $F(-2.45\,\,\text{GeV}^{2})=0.167\pm 0.010_{-0.007}^{+0.013}$. It is important to emphasize that our phenomenological input is complementary to that used in the determinations quoted by PDG. Finally, the condition $F(0)=1$ imposed by gauge invariance is implemented in an exact way.

Calculation of $\langle r_{\pi}^{2}\rangle$ and its uncertainty.— We assume first that only one spacelike value $F(t_{a})$ at a point $t_{a}<0$, and one value of the modulus $|F(t_{b})|$ at a point $t_{b}$ in the specified range of the elastic region $(4m_{\pi}^{2},t_{\text{in}})$ are known. From the input described above, it is possible to derive an exact range for the first derivative of $F(t)$ at $t=0$, equal by definition to $\langle r_{\pi}^{2}\rangle/6$. We do not give the proof, which can be found in [8, 9], but simply quote the result. In order to proceed, we must introduce some notation. First make a change of variable from $t$ to $z=\tilde{z}(t)$, where

[TABLE]

and define a new function $g(z)$ by

[TABLE]

Here $\tilde{t}(z)$ is the inverse function of $\tilde{z}(t)$, ${\cal{O}}(t)$ is the Omnès function

[TABLE]

where $\delta(t)$ is equal to the $\pi\pi$ $P$-wave phase-shift $\delta_{1}^{1}(t)$ for $t\leq t_{\text{in}}$ and is an arbitrary smooth (Lipschitz continuous) function above $t_{\text{in}}$, and $w(z)$ and $\omega(z)$ are two auxiliary functions, defined as

[TABLE]

[TABLE]

Then one can prove the following exact inequality, expressed as positivity of the determinant:

[TABLE]

where $z_{a}=\tilde{z}(t_{a})$, $z_{b}=\tilde{z}(t_{b})$ and $\bar{g}(z)=g(z)-g(0)-g^{\prime}(0)z$.

As shown in [9], the dependence of the Omnès function ${\cal{O}}(t)$ on the arbitrary phase $\delta(t)$ above $t_{\text{in}}$ is exactly compensated in the product (4) by the corresponding dependence of $\omega(z)$. Therefore, except $\langle r_{\pi}^{2}\rangle$, all the quantities entering the determinant (8) are calculable and depend only on the given input. From (4) it follows that $\langle r_{\pi}^{2}\rangle$ enters only the expression of the derivative $g^{\prime}(0)$, which can be written as $g^{\prime}(0)=\xi+\eta\,\langle r_{\pi}^{2}\rangle$, where $\xi$ and $\eta$ are real numbers. This implies that the determinant is a quadratic concave expression of $\langle r_{\pi}^{2}\rangle$ and the positivity condition (8) can be written as

[TABLE]

This inequality leads to a definite allowed range for $\langle r_{\pi}^{2}\rangle$ if $B^{2}-AC\geq 0$ and has no solution if $B^{2}-AC<0$. The latter case occurs when the phenomenological input adopted is inconsistent with analyticity.

The inequality (8) involves only one spacelike value $F(t_{a})$ and one timelike modulus $|F(t_{b})|$, which have been assumed to be known. The optimization formalism applied here allows actually the simultaneous inclusion of several spacelike and timelike points. The general case, treated in [9], is expressed by a more general determinant, but leads to a condition of the same form (9) for the charge radius. However, the simultaneous inclusion of many points leads to a system over-constrained by analyticity, which is difficult to handle numerically. Therefore, the most suitable approach for phenomenological applications is that adopted here, in which only one spacelike and one timelike constraint are imposed at the same time, the results from different inputs being subsequently combined by a suitable averaging procedure.

From the inequality (9), an allowed range for $\langle r_{\pi}^{2}\rangle$ can be obtained for every set of inputs. However, except the condition $F(0)=1$, which is exact, the input quantities are known only with some errors. One of the key aspect of our calculation is a proper statistical treatment of these errors. This is achieved by randomly sampling each of the inputs with specific distributions: the phase of $F(t)$, which is the result of a theoretical calculation, is assumed to be uniformly distributed, while for the spacelike and the timelike data, which are known from experimental measurements, we adopt Gaussian distribution with the measured central value as mean and the quoted error (the biggest error for spacelike data where the errors are asymmetric) as standard deviation.

For each set of values of the input statistical sample, if they are compatible, we calculate upper and lower value on $\langle r_{\pi}^{2}\rangle$ from (9). Since all the values between these extreme points are equally allowed, we uniformly generate values of $\langle r_{\pi}^{2}\rangle$ in between the bounds. For convenience, the minimal separation between the generated points was set at $10^{-3}\,\text{fm}^{2}$ and for intervals smaller than this limit no intermediate points were created. In this way, for each input from a spacelike energy and one timelike point in the region $(0.65-0.71)\,\text{GeV}$, we obtain a large sample of values of $\langle r_{\pi}^{2}\rangle$. The results were proved to be stable against the variation of the size of the random sample and the minimal separation mentioned above.

In Fig. 1 we show for illustration two distributions of $\langle r_{\pi}^{2}\rangle$, obtained using the Bern phase and one timelike modulus taken from the BABAR experiment [16]. Similar results have been obtained with the Madrid phase and other data on modulus. The left panel shows the distribution obtained without the inclusion of a spacelike datum (this case is obtained from Eq. (8) by removing the second row and column of the determinant). The histogram is rather flat and far from a normal distribution, which means that the input does not allow the extraction of a precise value of $\langle r_{\pi}^{2}\rangle$ and its error. Adding the input from the spacelike point $t=-1.60\,\text{GeV}^{2}$, we have obtained the distribution in the right panel, shown actually to be very close to a Gaussian. The explanation is that, adding the spacelike information, the input is much more constrained and many points in the Monte Carlo generated input do not satisfy the compatibility condition discussed below Eq. (9). The consistent values of the input lead to a normal distribution, which allows the extraction of the mean value and the standard deviation for the parameter $\langle r^{2}_{\pi}\rangle$.

In Fig. 2 we show the 68.3% confidence limit (CL) intervals obtained in this way, for all the input data on modulus in the chosen range $(0.65-0.71)\,\text{GeV}$, for all experiments. The Bern phase and the spacelike input at $t=-1.60\,\text{GeV}^{2}$ have been used in this figure. Similar results are obtained with the Madrid phase, and also using the spacelike datum at $t=-2.45\,\text{GeV}^{2}$.

We take then the average of the results obtained with input from various measurements. Since the degrees of correlations between the measurements at different energies are expected to vary from one experiment to another, we perform first the average of the values obtained with input from each experiment. The correlations are actually not known, therefore we apply the prescription proposed in [32] and adopted by PDG [1], where the effective correlation is determined from data themselves. Denoting by $v_{i}$ and $\sigma_{i}$ the central values and the standard deviations obtained from $n$ measurements, we define a generic covariance matrix $C(f)$ by $(C(f))_{ii}=\sigma_{i}^{2}$ and $(C(f))_{ij}=f\sigma_{i}\sigma_{j}$ for $i\neq j$, where $f\in(0,1)$ is an unknown global correlation. Then the most robust average $\bar{v}$ is defined as standard mean of $v_{i}$ weighted with the normalized weigths $\sigma_{i}^{-2}/\sum_{j}\sigma_{j}^{-2}$, and the standard deviation $\bar{\sigma}$ is given by

[TABLE]

where $f$ is the solution of the equation $\chi^{2}(f)/(n-1)=1$, for the standard definition of $\chi^{2}(f)$ [32].

We have applied the above prescription to the values obtained from each $e^{+}e^{-}$-annihilation and $\tau$-decay experiment, separately for the Bern and the Madrid phase. A conservative way of treating the information from spacelike $t$ has been adopted, by taking, for each timelike input, the simple average of the central values and standard deviations obtained with each of the two considered data. It turned out that in all cases the ratio $\chi^{2}(0)/(n-1)$ was less than 1 and increased for a positive correlation, reaching unity for $f$ in the range $0.68-0.96$.

Results.— In Table 1 we show the means and standard deviations obtained with the averaging prescription described above for the different timelike experiments. The results are mutually consistent among them, reflecting the consistency of the input data on modulus in the selected energy range. The spread of the errors, ranging from the smallest (BABAR) to the largest (OPAL), is actually not very large, which can be explained by the common information that entered the predictions. We can therefore combine all the results, without the risk of a bias from a single experiment. Since the correlations between these values are difficult to estimate, we have applied again the averaging procedure described above [32], where an empirical overall correlation is extracted from the data. The correlation for all the ten values in Table 1 was found to be 0.82 and 0.75 for the Bern and the Madrid phase, respectively, very close in both cases to the value $f$ where $\bar{\sigma}$ defined in (10) reached its maximum. For a conservative estimate of $\bar{\sigma}$ we have adopted actually the maximum value of (10) for $f$ in the range (0,1). This led to the predictions $0.4317\pm 0.0044\,\text{fm}^{2}$ and $0.4323\pm 0.0039\,\text{fm}^{2}$ for $\langle r_{\pi}^{2}\rangle$, for the Bern and the Madrid phase respectively. Once again adopting a conservative treatment, we have combined the two determinations by taking their simple average $\langle r_{\pi}^{2}\rangle=(0.4320\pm 0.0041)\,\text{fm}^{2}$, from which we obtained our final result

[TABLE]

The separate predictions obtained from $e^{+}e^{-}$ and $\tau$-decay data are $(0.657\pm 0.003)\,\text{fm}$ and $(0.658\pm 0.004)\,\text{fm}$, respectively, showing that (11) is dominated by the more precise $e^{+}e^{-}$ data.

The central value (11) is lower than the PDG average (1), a feature that seems to be common to the determinations based on analyticity [6, 7]. We note however that (11) is consistent with the determinations based only on elastic $e\pi$ scattering data, without including the more model-dependent determinations from electroproduction data at low spacelike $t$ [1]. If not settled by experimental considerations, this difference is a challenge for precise lattice calculations, expected in the coming years.

The present work proves the strength of the general principles of analyticity and unitarity, enforced by suitable mathematical techniques and supplemented by Monte Carlo simulations for error assessment. Our prediction (11) is based on input complementary to that entering the PDG average (1) and has an uncertainty smaller by a factor of almost 2.7, making it it the most precise model-independent determination to date.