The size of the pion from full lattice QCD with physical $u$, $d$, $s$ and $c$ quarks

We present the first calculation of the electromagnetic form factor of the $\pi$ meson at physical light quark masses. We use configurations generated by the MILC collaboration including the effect of $u$, $d$, $s$ and $c$ sea quarks with the Highly Improved Staggered Quark formalism. We work at three values of the lattice spacing on large volumes and with $u$/$d$ quark masses going down to the physical value. We study scalar and vector form factors for a range in space-like $q^2$ from 0.0 to -0.1 $\mathrm{GeV}^2$ and from their shape we extract mean square radii. Our vector form factor agrees well with experiment and we find $\langle r^2 \rangle_V = 0.403(18)(6) \,\mathrm{fm}^2$. For the scalar form factor we include quark-line disconnected contributions which have a significant impact on the radius. We give the first results for SU(3) flavour-singlet and octet scalar mean square radii, obtaining: $\langle r^2 \rangle_S^{\mathrm{singlet}} = 0.506(38)(53) \mathrm{fm}^2$ and $\langle r^2 \rangle_S^{\mathrm{octet}} = 0.431(38)(46) \mathrm{fm}^2$. We discuss the comparison with expectations from chiral perturbation theory.


I. INTRODUCTION
The electromagnetic form factor of the charged π meson parameterises the deviations from the behaviour of a point-like particle when struck by a photon. These deviations result from the internal structure of the π i.e. its quark constituents and their strong interaction. The form factor is calculable in QCD but a fully nonperturbative treatment is necessary at the small (negative) values of 4-momentum transfer, q 2 , covered by direct modelindependent experimental determination of the vector form factor [1] from π − e scattering. The experimental error is 1-1.5% in the region up to |q 2 | = 0.1GeV 2 and so a lattice QCD calculation of the form factor there can provide a stringent test of QCD. This is complementary to tests of QCD through calculation of meson masses and of decay constants that parameterise meson annihilation, for example to a W boson [2,3].
In the nonrelativistic limit, where q 2 ≈ −( q) 2 , the vector form factor, f + (q 2 ), can be viewed as the Fourier transform of the electric charge distribution. The mean squared radius obtained by integrating over this distribution is then given by This is adopted more generally as a definition of r 2 , since it is useful to have a single number with which to characterise the shape of the form factor. We will use it here to compare the 'size' of the π derived from our form factor with that obtained from experiment. * jonna.koponen@glasgow.ac.uk † christine.davies@glasgow.ac.uk ‡ http://www.physics.gla.ac.uk/HPQCD The calculation of r 2 from lattice QCD is complicated by the fact that the result is very sensitive to the mass of the π. The mean square radius diverges as m 2 π → 0 when the π meson cloud surrounding the π becomes of infinite range [4]. It has been numerically too expensive until recently to include u/d quarks with their physically very light masses in lattice QCD calculations. Results have instead had to be extrapolated to the physical point from heavier masses using chiral perturbation theory. The lightest π meson mass used in earlier calculations of the electromagnetic form factor has been in the range of 250-400 MeV, i.e. approximately twice the physical value or more. Results range from multiple values of the lattice spacing including the effect of u and d sea quarks [5][6][7] to those with only a single lattice spacing including the effect of u, d sea quarks [8] or, more realistically, u, d and s sea quarks [9][10][11]. See also [12] for a calculation in the ε regime. A mean square radius can similarly be defined for the scalar form factor. Earlier results, again for relatively heavy values of the π meson mass, have been obtained with u and d sea quarks (only) in [8,13,14].
Here we give results for both vector and scalar form factors for π mesons made of physical u/d quarks and including a fully realistic quark content in the sea, with physical u, d, s and c quarks. We also work with three different values of the lattice spacing. This enables good control of systematic errors both from m π and from discretisation. Our lattices have large volumes with a minimum spatial size of 4.8 fm.
The vector form factor is accessible in experiment and, as we shall see, our results can be directly compared to the experimental data with no extrapolation.
The scalar form factor cannot be obtained directly by experiment but information on it can be extracted by applying chiral perturbation theory to the π decay constant and to π − π scattering [4,15,16]. An additional ingredient in the lattice QCD calculation in this case is the need  [20,21]. The lattice spacing, a, is determined using the w0 parameter [23], and has a correlated 0.5% uncertainty from the physical value of w0, fixed using fπ [3]. Set 1 will be referred to "very coarse", 2 as "coarse" and 3 as "fine". Columns  to include quark-line disconnected contributions. The expectation from chiral perturbation theory [17] is for the disconnected contribution to the form factor at q 2 = 0 to be small but for the impact on the radius as defined in eq. (1) to be substantial. Our results are very much in line with expectations from chiral perturbation theory and we are able to distinguish disconnected contributions coming from u/d and s quark loops.
Section II describes how the lattice calculation is done and gives details of the results. Our results are compared to experiment, to chiral perturbation theory expectations, and to other lattice calculations in Section III and Section IV gives our conclusions, looking forward to improved calculations in future.

II. LATTICE CALCULATION
For the lattice QCD calculation we use the Highly Improved Staggered Quark (HISQ) action [18], which has been demonstrated to have very small discretisation errors [3,19]. We use gluon field configurations generated by the MILC collaboration [20,21] that include u, d, s and c sea quarks using the HISQ action along with a fully O(α s a 2 ) improved gluon action [22]. The ensembles that we use here have light quark masses m u = m d = m l with m l and hence m π close to its physical value. The parameters of the ensembles are given in Table I.
On these configurations we generate HISQ light quark propagators with the same mass as that of the sea light quarks. We use a local random wall source [3] and 4 time sources per configuration for high statistics. The propagators are combined into π meson correlation functions (2-point correlators) that create a π meson at time 0 and destroy it at time t and correlation functions that allow for interaction with a current J at an intermediate time, t, between a π meson source at 0 and sink at T (3-point correlators). These are illustrated in Fig. 1. Results at all t (t) values are obtained for 2(3)-point functions and we also use three values for T in the 3-point functions, so that our fits can map out fully the t and T dependence for improved accuracy. When J is a vector current we need to consider only one 3-point diagram for the flavour non-singlet π. This is shown as the central diagram of Fig. 1 in which the current J is inserted into one of the legs of the 2-point function. We simply multiply by 2 to allow for its insertion into the other leg. The 'disconnected diagram' which is the product of a π 2-point function and a closed quark loop coupled to J is shown as the lower diagram of Fig. 1. This vanishes for vector J in the ensemble average because it is odd under charge-conjugation [24]. For scalar J this diagram needs to be included and different combinations of flavours of quarks in the closed quark loop give rise to different form factors. The π mesons in our correlators are the Goldstone mesons whose mass vanishes with m l . We ensure this by using the local γ 5 operator at source and sink. In staggered quark parlance this is the γ 5 ⊗γ 5 operator. For J we use a symmetric 1-link point-split spatial vector current, V i , or a local scalar current, S. A gluon field is included in the vector current to make it gauge-covariant. Both of these are 'tasteless' staggered quark operators (γ i ⊗ 1 and 1 ⊗ 1) and so can be used in a 3-point function with the Goldstone meson at source and sink.
Results for the ratio of the 3-point correlator containing a vector current to the product of appropriate 2-point correlators for the pion on fine ensemble set 3. The ratio for the 3 different values of T are plotted as a function of t with a π of momentum zero on the left, and momentum ap = 0.0363 on the right. Note that this figure is to illustrate the quality of our results; we do not use this ratio to extract ground-state parameters. Instead we perform a simultaneous fit to multiple exponentials for both the 2-point and 3-point correlators as described in the text.
We work with several π meson spatial momenta by generating light quark propagators with a phase included on the spatial gluon links. This is equivalent to introducing a phase into the boundary condition on the field [25], which gives a momentum to the quark. This is referred to as using 'twisted boundary conditions'. As illustrated in the central diagram of Fig. 1 we choose the spectator quark in the 3-point function to have zero momentum and give momenta p 1 and p 2 to the quarks that interact with the current. Both momenta are chosen to be in the (1, 1, 1) direction. By using various values of p 1 and p 2 we can obtain 3-point functions at several different small values of squared 4-momentum transfer, q 2 .

A. Vector Form Factor
For the vector current case we have a set of 2-point and 3-point quark-line connected correlators at various values of p 1 and p 2 on each ensemble. The quality of our results is illustrated for one ensemble and set of momenta in Figure 2. The 2-point and 3-point correlators are all fit simultaneously using Bayesian methods [26] that allow us to include the effect of excited states, both 'radial' excitations and, because we are using staggered quarks, opposite parity mesons that give oscillating terms. Since the oscillating terms are absent for zero-momentum π mesons they are small here, but we nevertheless include them in our fits. Having 3-point correlators from multiple T values is also important in taking account of excited states. Fitting multiple momenta simultaneously allows us to take account of correlations between the correlators.
The fit form for the 2-point function with source at 0 and sink at t and spatial momentum, p, is: Opposite parity terms (o.p.t.) are similar to the terms given explicitly above but with factors of (−1) t /a . For the 3-point function [27,28]: Prior values and widths are taken as: ground-state energy, 10% width; splitting between ground-state and excited energies, 650 MeV with 50% width; splitting between ground-state and lowest oscillating state, 500 MeV with 50% width; amplitudes, 0.01(1.0) for normal states and 0.01(0.5) for oscillating states; matrix elements, 0.01(1.0) for vector currents. We take the result from a 6 exponential fit (with 6 oscillating exponentials) to obtain the vector form factor. The ground-state parameters are given by i = j = 0 in eqs. (2) and (3) and are our key results. By matching to a continuum correlator with a relativistic normalisation of states and allowing for a renormalisation of the lattice current, we see that the matrix elements between the ground state mesons that we want to determine are given by: The matrix element is related to the form factor for the vector current via: The vector matrix element can be normalised using the fact that f + (0) = 1 for a conserved current (inserted in either the quark or the antiquark legs in Fig. 1), and we can therefore determine Z by demanding that condition for our current. f + (0) is determined at q 2 = 0 for spatial V i by setting p 1 = p 2 = 0. Table II gives results for the π energies and 2-point amplitudes as a function of momentum. A good test of discretisation errors is to determine the speed of light, c 2 from (E 2 − m 2 )/p 2 . From Table II we see that c 2 deviates from 1 at most by 2(1)% at the largest momenta. Another test is to compare the scaling of amplitudes to the expected 1/ √ E behaviour for a pseudoscalar. Again we see good agreement, with deviations at most 3(1.5)%. Results for π energies in lattice units at the different spatial momenta used on each set, as well as the corresponding amplitudes from the 2-point functions. The values given here come from the simultaneous fit of 2-point correlators with 3-point correlators containing a vector current. Results for 2-point parameters from the fit of 2-point correlators with 3-point correlators including a scalar current are in agreement. π results at zero momentum agree with those in [3,29], but are not the same because the fits used here also include 3-point functions. Lower table: Results for unrenormalised form factors at q 2 values corresponding to different combinations of π momenta (from the upper table) at source and sink. The results at q 2 = 0 come from using the lowest non-zero spatial momentum for both p1 and p2. The scalar form factor results given are for the connected 3-point function only. Table II also gives the raw (unrenormalised) form factors for various q 2 values obtained from different combinations of momenta (in positive and negative) directions at source and sink. The statistical errors on the form factors are 0.5-3%. By dividing the values at non-zero q 2 by the value at q 2 = 0 we obtain normalised values for f + . f + is plotted against q 2 in Fig. 3 for all three sets along with the results from experiment [1]. The agreement with experiment is good, reflecting the fact that our results correspond to physical π masses.

Results
In fitting a functional form in q 2 to our results to extract a mean squared radius, we use the same form as that used for the experimental results [1], but including allowance for finite lattice spacing effects. We also fit over a similar range of q 2 values. We use: Here c a n and b a n allow for discretisation effects in the normalisation of f + and in r 2 respectively. We take Λ = 500 MeV and allow priors on the b and c fit parameters of 0.0(1.0) for b a 4 and c a 4 and 0.0(0.3) on b a 2 and c a 2 (since tree-level a 2 errors are absent in this calculation). We allow a prior width on the physical result for the mean squared radius, r 2 V of 25%. The term with coefficient b sea allows for mistuning of sea quark masses. From chiral perturbation theory a term linear in the quark masses is expected, and it is convenient to take this term as a ratio to another quark mass so that factors of the quark mass renormalisation cancel. The factor of 10 multiplying m s,phys gives a value close to the chiral scale, Λ χ . The mistuning of the sea masses is defined as ) and values of δm sea /m s,phys values for these ensembles are tabulated in [30]. The values are all less than 0.05, but not zero because of mistuning of the sea s quark mass.
The final logarithmic term in eq. (7) comes from chiral perturbation theory [4] and is the source of the divergence in the radius as m π → 0. We use it, rescaling the argument of the logarithm so that it vanishes at the physical pion mass, to make small adjustments for the fact that our u/d quark masses are not exactly at their physical values (in fact they are slightly too low). Λ χ = 1.16 GeV. Because we are very close to the physical light quark mass point we do not need to include further terms in a chiral perturbation theory expansion since they will be negligible.
We apply the functional form of eq. (6) and (7) to our result taking account of the correlations between results at different values of q 2 obtained on a given ensemble. The fit has χ 2 /dof = 0.9 and gives the physical result for the electric charge radius of the π of r 2 V = 10.35 (46) We can also use the final logarithmic term in eq. (7) to estimate the impact of isospin and electromagnetic effects by varying the value of m π,phys used there. The physical value of m π corresponding to our lattice world in which u and d quark masses are equal and there is no electromagnetism is m π 0 = 0.135 GeV [31], and we use this for our central value above. The experimental results correspond to m π + = 0.139 GeV and we substitute that for the physical value in the logarithm to assess the uncertainty from the fact that the real world has different u/d quark masses and the quarks have electric charge. This gives an estimate for the systematic uncertainty from isospin/electromagnetism of 0.5%.
We must also include a systematic uncertainty from working on lattices with finite spatial volume, albeit large. Finite volume effects are small on these lattices for the π mass and decay constant [3] and effects of similar size are expected in the form factor at fixed q 2 . Because the mean squared radius is defined from the small difference in values for the form factor as q 2 moves away from zero (where the form factor is defined to be 1), a small effect on the form factor at non-zero q 2 can become a significant effect on the radius. These effects can be estimated from chiral perturbation theory. Continuum chiral perturbation theory is a good guide here and we do not need staggered chiral perturbation theory because, as shown in [32], staggered quark taste-effects which might be expected to affect π masses appearing in chiral loops in fact tend to cancel against associated hairpin diagrams. It turns out that this cancellation happens for a wide range of quantities (including decay constants and form factors) for a specific value of the hairpin coefficients that seems to be close to the value obtained in practice. We therefore use continuum analyses and specifically results from analyses that are relevant to our use of twisted boundary conditions [33,34] because this modifies the expected finite-volume dependence. From [33] the relative finite volume effect in the vector squared-radius varies in the range 1-1.5% for lattice sizes that we use in the range 4.8 fm to 5.8 fm for physical π masses. Note that the direction of the finite-volume effect is such that the radius would be larger in the infinite volume limit. We do not make a correction for this but include an uncertainty of 1.5% for finite volume effects.
Our error budget for r 2 V is given in Table III. Adding the systematic uncertainties in quadrature as the second uncertainty gives our result: to be compared to 0.431(10) fm 2 from the experimental results of [1] using the same fit form. The Particle Data Group [35] give a mean square radius from averaging over several experimental results of 0.452(11) fm 2 .
B. Scalar Form Factor

Results for the connected contribution
We begin by discussing our results for the connected contribution to the scalar form factor of the π. This is the result calculated from 3-point functions of the form sketched in the central diagram of Fig. 1 in which the scalar current is composed of the light quarks which are the valence quarks of the π. Although this form factor does not correspond to a physically realisable process (even if we had a particle with which to produce a scalar current) it is nevertheless possible and useful to compare different lattice QCD calculations for it. Different formalisms within lattice QCD should give the same results in the continuum and chiral limits for the mean square radius from the connected scalar form factor. A key issue, to be discussed further below, is then how big the additional contribution is from quark-line disconnected diagrams.
The calculation for the connected scalar form factor proceeds in an identical way to that of the vector form factor discussed in Section II A. We calculate the 3-point function given as the central figure of Fig. 1 with a scalar current made from light quarks inserted as J. We use the same light quark propagators and 2-point functions as for the vector case. The quality of our results is illustrated for one ensemble and set of momenta in Fig. 4 (we use the same ensemble and set of momenta as in Fig. 2).
We fit the 2-point and 3-point correlators simultaneously (but in a separate fit from the vector case) as a function of t, t and T as given in eqs. (2) and (3). The priors are taken to be the same as in the vector case except that the prior width on the scalar matrix element is taken to be much larger, reflecting expectations on its value given below. We take the prior width on the scalar The ground-state matrix element for the scalar current is related to the parameter J 0,0 extracted from our fits as in eq. (4). In turn the matrix element is related to the form factor that we wish to extract by where A is a normalisation factor. Our scalar current made from HISQ quarks is absolutely normalised [36]. If we had included disconnected diagrams associated with the scalar current we would be able to write, from the Feynman-Hellmann theorem, with f 0 (0) = 1, A = (∂m 2 π /∂m l )/2 and we take the same m l value for valence and sea l quarks. The factor of 2 on the right-hand side comes from the fact that we are inserting a scalar current in only one propagator to make the 3-point correlator and the π has two valence light quarks. For the connected correlator we expect instead that factor A in eq. (9) should be equal to half the derivative of the squared π mass with respect to its valence quark mass. This can be tested approximately using π and K masses on these ensembles in [3]. Comparing a(m 2 K − m 2 π )/(m s − m l ) (i.e. an approx derivative for a pseudoscalar meson mass) to the result for the unrenormalised f conn 0 (0)/Z S (=A) from Table II shows agreement within 10%, confirming that A does have the expected value. Since here we are chiefly concerned with the shape of the form factor, we simply treat the scalar current as requiring a Z factor, Z S , and determine this from the requirement that also f 0 (0) = 1.
To determine the mean squared radius associated with the connected scalar form factor we take the same fit as for the vector case, eqs. (6) and (7), except that the coefficient of the chiral logarithm is now 6. This coefficient applies to the complete scalar form factor but we use it here to estimate conservatively the impact of changing the π mass close to the physical point and therefore the uncertainty. We use the same priors on the coefficients as in the vector fit, except that we increase the prior on the physical result for the mean squared radius, since its value is less well-known.
Our fit has a χ 2 /dof of 1.1 and gives a physical result, r 2 conn S of 8.97(45) GeV −2 or 0.349(18) fm 2 . The systematic errors are somewhat larger in the scalar case because of the larger coefficient of the chiral logarithm. Using this we estimate the systematic uncertainty from isospin/electromagnetism at 3%. The larger coefficient for the chiral logarithm also carries with it the implication of larger finite volume effects, potentially by a factor of 6, giving a systematic uncertainty of 10% on our ensembles from this source, allowing for the fact that the mean square radius is slightly smaller than for the vector case.
Adding systematic errors in quadrature our final result for the mean squared radius from the connected scalar form factor is r 2 (π) S,conn = 0.349(18)(36)fm 2 Our error budget is given in Table III. This radius has a central value that is only slightly smaller than the vector form factor radius (eq. (8)). This is illustrated in Fig. 5  in which we compare the lattice results and the fit results for the vector and connected scalar cases.

Results including the disconnected contribution
For the full scalar form factor we need to include the quark-line disconnected contribution from the lower diagram of Fig. 1. We can then define flavour-singlet and flavour-octet scalar currents: S octet = 2ll − 2ss.
Scalar form factors for these two currents are then determined by combining the connected scalar form factor of Section II B 1 with disconnected contributions in appropriate combinations from quark loops made from l quarks or s quarks. For the q 2 = 0 case it is particularly simple to calculate the disconnected contributions. Indeed, for the ss scalar current this is the π meson equivalent of the 'strangeness in the nucleon' calculation on which there has been a great deal work in lattice QCD (see, for example, [37]). The disconnected contribution for current qq is π|S qq |π disc = π(p)|qq|π(p) − π(p)|π(p) qq . (13) The first term is the ensemble average of a π meson 2point function with source at time 0 and sink at time T with a scalar current (condensate) insertion summed over the spatial points making up the timeslice at t. The second term in eq. (13) subtracts the product of the vacuum expectation values of the π meson correlator and condensate.
With HISQ quarks a convenient way to represent the qq condensate is as a sum over a pseudoscalar meson correlator with valence quarks q [29]. We use an identity [29,38] that relates the quark propagator for staggered quarks on a given gluon field configuration to a product of quark propagators summed over lattice sites: Here 0 and n are arbitrary lattice sites, Tr is a color trace and am q is the quark mass in lattice units used in the propagator. The left-hand side of eq. (14) is the negative of the quark loop from site 0 to site 0 needed for the disconnected piece of the scalar current and the the righthand side is the Goldstone pseudoscalar meson correlator at zero spatial momentum for a quark-antiquark pair of mass am q multiplied by that mass. Since our Goldstone meson correlators here use a random-wall source a sum over a timeslice for lattice site 0 for the quark loop is done implicitly. The quantities required to calculate the disconnected contribution for the ss current to the scalar form factor at q 2 = 0 are then simply π meson correlators and those for the pseudoscalar ss meson known as the η s . The η s does not correspond to a physical particle but its correlators are nevertheless usefully studied in lattice QCD [39] and so have been calculated previously [3]. To make the 3-point function needed we take a π meson correlation function with source at timeslice 0 and sink at time-slice T and a set of η s meson correlators with sources at timeslices denoted by t. For this calculation we use correlators made for the determination of π, K and η s masses and decay constants in [3]. These have zero spatial momentum and 16 time-sources, so that t comes in steps of 4 timeslices on coarse set 2, which is the set that we will focus on. The ss current loop at time t is obtained by summing over all end points for an η s correlator starting at time-source t. For the disconnected contribution we multiply this by the π correlator with time-source 0 and sink T , averaging over all time-sources on a configuration that give the same set of relative time separations. The 3-point function that yields π|S ss |π of eq. (13) at q 2 = 0 is thus given, averaging over gluon field configurations, by where p 1 = p 2 = p = 0. For the scalar current made of light quarks an equivalent expression holds, using two π meson correlators with offset time-sources. Figure 6 shows results for a scalar current made of light quarks or strange quarks. The quantity plotted is the ratio of the 3-point correlator generated from the equivalent of eq. (15) divided by the 2-point correlators for the π meson at zero momentum whose sources are at T lattice spacings apart. We take T = 32 but have checked that results are very similar for other values of T , such as T = 28. Because we are using point sources for our π meson correlators we do not have a large plateau region. A longer plateau is obtained using smeared sources [40]. However, by using a combined fit of the 3-point and 2point correlators we can allow for systematic uncertainties from excited states and we obtain a good fit. The red and blue hashed bands show the fit results for the (unrenormalised) ground-state matrix element of the scalar current made of light and strange quarks, respectively, divided by twice the π meson mass. The results for the ground-state matrix elements for both the l and s scalar currents are given in Table IV. At q 2 = 0 we can compare them to the connected contribution given in Table II. We see that the disconnected contributions are very much smaller, each around 1% of the connected contribution. This is to be expected based on the Feynman-Hellmann theorem which would relate the disconnected contributions to the derivative of the π meson mass with respect to the sea s or l quark mass [37], in a similar way to that discussed in Section II B 1 for the connected contribution. Our results for the s scalar current indicate reasonable agreement (within a factor of two) with the π mass dependence on the s sea quark mass (keeping all other parameters fixed) obtained at heavierthan-physical π masses by the MILC collaboration [40].
To obtain results for the disconnected contribution to the scalar form factor at non-zero values of q 2 we need to project onto non-zero lattice spatial momenta, 2π/L s (n x , n y , n z ), at T and t in the correlators used in eq. (15). This extends eq. (14) to the non-zero momentum case using translation invariance. The statistical errors grow as spatial momentum is introduced so we restrict ourselves to the smallest non-zero lattice momenta with (n x , n y , n z ) equal to (1, 0, 0) and (1, 1, 0) and permutations thereof, including −1 as well as +1. We work only on coarse set 2.
We obtain the ground-state matrix element for these contributions using a combined 3-point and 2-point fit as before. For this we need new 2-point π meson correlators at these spatial momentum values and we obtain these on a subset of 600 configurations. Figure 7 shows the disconnected contribution to the scalar form factor for l and s scalar currents, tabulated in Table IV. It is clear from this plot that when disconnected contributions are included with a positive sign they will increase the slope of the form factor and therefore the  (12)). These are plotted as a function of q 2 for coarse lattices, set 2.
mean square radius. Thus the flavour singlet scalar radius will be larger than the radius from the connected diagram only. This is less clear for the flavour octet case but the magnitude of the strange disconnected contribution is smaller than that of the light disconnected contribution so we might expect a net positive effect. Figure 8 collects the results into singlet and octet flavour combinations. Now it is clear that both singlet and octet radii will be larger than the radius from the connected diagram only.
To obtain the mean square radius for the singlet and octet scalar form factors we must combine the connected and disconnected contributions. In doing this we must be careful to insert appropriate factors of 2 for the ll pieces so that both the connected and disconnected contributions include uu + dd. Since we only have a calculation of the disconnected pieces on coarse set 2 we use a simple approach to determining the change in the mean square radius, using a linear approximation to the form factor over the small q 2 range (0 to -0.0315 GeV 2 ) covered by the disconnected results. This has the advantage of making clear how the disconnected contributions affect the result. They appear both in the value of the total form factor at q 2 = 0 which is used for the normalisation and they contribute to the slope of the form factor in q 2 . As discussed above, the effect on the form factor at q 2 = 0 is very small (1%) and the largest effect comes from the contribution to the slope. We have, comparing the form factor at q 2 to that at 0, The second term makes a large contribution to the mean square radius because the change in the disconnected contribution to the form factor over the range in q 2 (depending on the combination of flavours) is of the same size as that of the connected contribution included in r 2 conn . We find, for example, that the change in mean square radius is 50(20)% for the singlet combination. For the singlet and octet combinations we obtain: r 2 (π) S,octet = 0.431(38)(46) fm 2 .
Here the first error is statistical and comes from adding in the disconnected contribution. The second error is systematic from electromagnetic/isospin and finite volume effects as discussed in Section II B 1 for the connected scalar radius. The full error budget for the singlet/octet radius is given in Table III. For comparison with earlier work on configurations that include only u and d quarks in the sea we can construct a radius that corresponds to the form factor for a uu + dd scalar current. We find S,ud = 0.481(37)(50) fm 2 .
As eq. (16) makes clear, the results for the different scalar radii are correlated. The differences between them are significant since a lot of the uncertainty cancels. For example S,singlet − r 2 (π) S,octet = 0.075(20) fm 2 .
We find the ordering: S,singlet > r 2 (π) S,ud > r 2 (π) S,octet > r 2 (π) S,conn . (20) FIG. 9. A summary of lattice QCD results for the mean square electric charge radius of the π meson arranged by the number of quark flavours included in the sea. The top result is from this paper; those including u, d, and s quarks in the sea (n f = 2 + 1) are from [9][10][11] and those including only u and d quarks in the sea (n f = 2) are from [5][6][7][8]. Results that include only one value of the lattice spacing have dotted error bars. Experimental results are from [1,[41][42][43][44]. The hashed vertical line gives the average from the Particle Data Group [35]. Figure 9 compares the result obtained in this paper for the mean square of the pion electric charge radius to other lattice QCD calculations by RBC/UKQCD [11], PACS-CS [10], the Mainz group [7], QCDSF [5], ETMC [6] and JLQCD/TWQCD [8,9], and to experimental results [1,[41][42][43][44]. It should be noted that several of these calculations include results at only one value of the lattice spacing and error budgets are not complete in all cases. A recent calculation by B. Owen et. al. [46] used one lattice spacing and five different pion masses down to 156 MeV but no chiral or continuum extrapolation is given so the results are not included in the figure.

III. DISCUSSION
The calculation presented in this paper is the first one that has been done at the physical pion mass -other lattice QCD calculations have used heavier than physical pions. However, as Figure 9 shows, all lattice QCD results agree well after extrapolation to zero lattice spacing and physical pion mass. We see no difference between the lattice calculations using different sea quark content (u and d only, u, d and s, or u, d, s and c quarks in the sea) at this level of accuracy. In Figure 3 we compare the shape of the electromagnetic form factor from our calcu- The HPQCD Collaboration's results are from this paper: "connected" shows the mean square radius from the quark-line connected calculation only (eq. (11)); "singlet" and "octet" are full calculations including quark-line disconnected diagrams arranged in flavoursinglet or flavour-octet currents (see eqs. (12) and (17)) and uu + dd includes only u/d quarks in the scalar current (eq. (18)). The results including only u and d quarks in the sea (n f = 2) are from [8,14]. The hashed green vertical bands give the result expected from chiral perturbation theory for Fπ/F for n f = 2 and n f = 2 + 1 (for comparison with our n f = 2 + 1 + 1 results) from eq. (22). The phenomenological result from π−π scattering [4] is very similar to the n f = 2+1 green band. The hashed purple vertical band gives the chiral perturbation theory expectation for the scalar octet mean square radius (eq. (25)) [45].
lation to the result by NA7 [1], which is the most accurate one of the experimental results. The agreement is good, which shows also here when we compare our result for the mean square radius to the NA7 results. Our value is 2σ below the average of the experimental results [35].
In the case of the scalar radius, comparison with other lattice QCD results must be done with care. Two lattice QCD calculations have been done including u and d quarks in the sea (n f = 2). These are by the Mainz group [14], recently updating [13], and the JLQCD/TWQCD collaborations [8]. Both of these calculations include the quark-line disconnected diagrams but only for a uu + dd scalar current (consistent within an n f = 2 framework). Here we include u, d, s and c quarks in the sea and a scalar current that includes also ss contributions in two different overall flavour combinations. We neglect cc contributions to the current since we expect those contributions to be suppressed by powers of the c quark mass.
The pion scalar form factor is not directly accessible to experiment as there is no suitable low-energy probe. However, the scalar radius can be determined from the cross section for π −π scattering and from the pion decay constant by using chiral perturbation theory [4,45]. In SU(2) chiral perturbation theory the ratio of the physical pion decay constant to that in the chiral limit can be related to the pion scalar radius for the uu + dd scalar current by [4] where F π = f π / √ 2=92 MeV and we take m π = 135 MeV. F π /F values from SU(2) chiral perturbation theory analyses of lattice QCD calculations are collected in [47] and give averages: (67), r 2 (π) S,ud = 0.76(9)(4)fm 2 , n f = 2 (22) = 1.0624 (21), r 2 (π) S,ud = 0.61(3)(3)fm 2 , n f = 2 + 1 where we have included a second uncertainty of 5% for higher order corrections to the chiral perturbation theory formula. The results for n f = 2 and n f = 2 + 1 analyses are compatible within 2σ but do not have to agree. A phenomenological estimate for r 2 (π) S,ud based on π − π scattering gives 0.61(4) fm 2 [16], which agrees well with the n f = 2+1 result above. For our n f = 2+1+1 calculations we compare our value for r 2 (π) S,ud to the n f = 2 + 1 results, following the discussion of the compatibility of n f = 2 + 1 and n f = 2 + 1 + 1 chiral analyses in [47].
When ss components are included in the current we can form flavour octet and singlet combinations. The flavour octet combination is interesting because it can be estimated from f K /f π since no new low-energy constants appear [45]. The mean square radius is given by with Using m π = 135MeV, m K = 496MeV, m η = 548MeV and F K /F π = 1.1916 [3] gives where we take the estimate from [45] of the uncertainty from higher order terms. Note that we expect this mean square radius to be smaller than r 2 (π) S,ud because it involves the subtraction of twice the strange current quarkline disconnected contribution, and our results show this to have a positive impact on the mean square radius.
The difference between singlet and octet mean square radii comes from the q 2 dependence of the matrix element of the ss piece of the scalar current. It is denoted δr 2 S in [45] and estimated in chiral perturbation theory as: improved Wilson HISQ FIG. 11. A comparison of the ratio of quark-line disconnected to connected contributions for the uu + dd current to the pion scalar form factor at q 2 = 0 plotted against the square of the lattice spacing. Results denoted by the blue plus symbol for an improved Wilson action are taken from [14]. The red triangles are from the results presented here for the HISQ action.
Using L r 4 renormalised at m η from [3] gives δr 2 S = 0.015 ± 0.1 fm 2 which, given its uncertainty, agrees with our result in eq. (19). Figure 10 compares our results for our various pion scalar mean square radii to those obtained on n f = 2 gluon configurations and to the expected values from chiral perturbation theory given above. There is reasonable agreement (within 2σ) between all the lattice QCD results for r 2 (π) S,ud and with the values expected from F π /F . Our result for the flavour octet radius is also in good agreement with the value in eq. (25). To illustrate how important the contributions from the disconnected diagrams are to the various scalar radii we also show the result from our calculation of the connected diagram only (eq (11)). As discussed in Section II B 2 the contributions from the disconnected diagrams to the form factor are small but the change in the slope and therefore in the radius is substantial. This feature of the scalar radius is also discussed in [13,14].
Our results are in both qualitative agreement and reasonable quantitative agreement with the picture expected from chiral perturbation theory [17]. There the disconnected contribution to the form factor is predicted to be very small at q 2 = 0, becoming negative at negative q 2 values so that the contribution to the radius is substantial (approximately equal to that of the connected contribution) and positive. We find the contribution to amount to an approximately 50% increase in the radius, rather than doubling, but with substantial uncertainty.
Although different lattice QCD formalisms will have a different normalisation for the scalar current, the ratio of disconnected to connected contributions to the form factor at a given value of q 2 should agree in the chiral and continuum limits since renormalisation factors will then cancel. Our results obtained here with the HISQ action seem to agree well with those from the overlap action given at one value of the lattice spacing in [8], judging this from their Figure 9. They do not agree well with those from an improved Wilson action given at three values of the lattice spacing in [14]. They have a very substantial disconnected contribution that also shows, as a ratio of the connected contribution, a very strong lattice spacing dependence. They work at heavier values of m π than we do but see little dependence on m π in this ratio (apart from one point that they suggest to treat as an outlier). In Figure 11 we show a comparison of their disconnected/connected ratio at q 2 = 0 to ours. For their points we have taken numbers from their Figure 2, using the values closest to m π = 250 MeV (so any variation with m π is not included in our plot, and it should be taken as an approximate representation of their results). For our results we give values from our coarse set 2 and from fine set 3. Although our results also show some lattice spacing dependence, it is hard to avoid the conclusion that the improved Wilson results are dominated by a lattice artefact. It is not clear whether the values from HISQ/domain wall and improved Wilson actions will agree in the continuum (and chiral) limits and this does need to be resolved.

IV. CONCLUSIONS
We have given the first lattice QCD results for the vector and scalar form factors of the π meson including u/d quarks with their physical masses.
Our results for the vector form factor as a function of q 2 agree well (needing no extrapolation) with the experimental values from π − e scattering (see Figure 3 and eq. (8)). This confirms the encouraging picture seen by earlier lattice QCD calculations (albeit with heavier u/d quark masses and a less realistic QCD vacuum) that lattice QCD does indeed reproduce the QCD effects that result in a finite electric charge radius for the π meson. It would clearly be possible to extend our results to higher values of q 2 with the aim of eventually matching on to expectations from QCD perturbation theory (see, for example, [48]). This would require finer lattices to avoid large discretisation errors from the discretisation of momentum and, for numerical efficiency, should probably be done with heavier-than-physical π mesons (even η s mesons).
The π scalar form factor is of less immediate phenomenological interest but can be stringent test of lowenergy expectations from QCD and from chiral perturbation theory (where it can be related to decay constants and π − π scattering). Here we have given the first results to include u, d and s quarks in the sea and in the scalar current. This allows us to define a number of dif-ferent radii for different flavour combinations. Calculation of the scalar form factor must include the effect of quark-line disconnected diagrams and we agree with earlier results that these have a substantial impact on the determination of the radii. An increase in the radius occurs where qq has a positive coefficient in the combination that appears in the scalar current and we find the magnitude of the contribution of ll to be larger than that of ss. We therefore have an ordering in value of the different radii that we give in eq. (20). Our values for the quark-line disconnected contribution to the form factor are, however, very small at q 2 = 0 in agreement with expectations from chiral perturbation theory and with earlier results using the overlap formalism [8].
Our value for the radius obtained using a u/d scalar current (eq. (18)) agrees within 2σ with expectations from chiral perturbation theory and earlier values from calculations including u and d quarks in the sea. We give the first results for the radius from the flavour octet and flavour singlet scalar currents (eq. (17)). The flavour octet mean square radius agrees well with expectations from chiral perturbation theory where it can be related to f K /f π and combinations of meson masses.
Our largest source of uncertainty in the scalar case is from finite-volume corrections since we are working with physically light π mesons, even if on lattices 5.8 fm across. To improve uncertainties on the scalar radius we would need to use larger volumes at the physical point, or include a dedicated study of finite-volume effects away from the physical point, and ensembles of gluon field configurations exist on which this could be done. Improved statistics on the quark-line disconnected contributions are also necessary to reduce the statistical uncertainty in their impact on the scalar mean square radius. The new techniques we have introduced here for handling the disconnected contributions in the pion scalar form factor make these improvements feasible in future results.