Model-independent constraint on the pion scalar form factor and light quark masses

We investigate the pion scalar form factor in the Meiman-Okubo framework, implementing the phase below the inelastic $K\bar K$ threshold, where it is known from the $\pi\pi$ scalar isoscalar phase shift $\delta_0^0$ by Watson theorem. State-of-the-art knowledge of the perturbative QCD expansion of the scalar correlator and the phase shift $\delta_0^0$ is used as input. No assumptions about the phase above the inelastic threshold or the possible zeros of the form factor in the complex plane are necessary. We obtain a model-independent constraint relating the sum of the light quark masses to the slope and the curvature of the pion scalar form factor at the origin. The recent lattice results for the light quark masses and the pion scalar radius are found to satisfy this constraint. We obtain also a strong correlation between the pion scalar radius and the curvature of the form factor, with rather high values predicted for the curvature.

We investigate the pion scalar form factor in the Meiman-Okubo framework, implementing the phase below the inelastic KK threshold, where it is known from the ππ scalar isoscalar phase shift δ 0 0 by Watson theorem. State-of-the-art knowledge of the perturbative QCD expansion of the scalar correlator and the phase shift δ 0 0 is used as input. No assumptions about the phase above the inelastic threshold or the possible zeros of the form factor in the complex plane are necessary. We obtain a model-independent constraint relating the sum of the light quark masses to the slope and the curvature of the pion scalar form factor at the origin. The recent lattice results for the light quark masses and the pion scalar radius are found to satisfy this constraint. We obtain also a strong correlation between the pion scalar radius and the curvature of the form factor, with rather high values predicted for the curvature.

I. INTRODUCTION
The pion scalar form factor Γ π (t) is defined by the matrix element π a (p)π b (p ) out|S(0)|0 = δ ab Γ π (t) , t = (p + p ) 2 , (1) of the scalar operator (2) where u, d are quark fields and m u , m d the quark current masses.
Since the Higgs boson is not light, the pion scalar form factor is not accessible to experiment. However, it is important for theory, reflecting crucial aspects of QCD at low energy. Its Taylor expansion at t = 0: Γ π (t) = Γ π (0) 1 + 1 6 r 2 π s t + c π s t 2 + . . . , convergent in a disk limited by the nearest branch point of Γ π (t) at t = 4m 2 π , has been investigated in chiral perturbation theory (χPT), where the pion scalar form factor has been calculated up to two loops [1][2][3].
The uncertainty in this relation might be somewhat underestimated, since Ref. [4] included only parts of the higher-order corrections. We shall discuss in the last section the impact of a larger uncertainty on the results derived in this paper. The quadratic scalar radius r 2 π s = 6Γ π (0)/Γ π (0) is connected to another important quantity of χPT, the effective chiral constantl 4 that determines the first nonleading contribution in the chiral expansion of the pion decay constant f π . It also contributes to the S-wave ππ scattering lengths a 0 0 and a 2 0 [5]. From general principles it is known that Γ π (t) is an analytic function of hermitian (real) type (i.e., it satisfies the Schwarz reflection relation Γ π (t) * = Γ π (t * )) in the t complex plane with a cut determined by unitarity for t ≥ 4m 2 π . Watson final-state theorem states that below the first inelastic threshold, which in practice is due to the KK channel, the phase of the form factor is equal to the phase-shift δ 0 0 (t) of the I = L = 0 partial-wave amplitude of ππ elastic scattering: There have been some discussions in the literature about the value of r 2 π s obtained in the frame of dispersion theory. While the treatments [5][6][7] based on Mushkhelishvili-Omnès equations give r 2 π s in the range (0.57 -0.65) fm 2 , the calculations [8,9] based on singlechannel Omnès formalism led to a higher prediction, r 2 π s = (0.75 ± 0.07) fm 2 . This discrepancy was discussed in [10], where it was shown that the single-channel treatment can be made consistent with the multi-channel one if one takes into account the fact that Watson theorem is valid modulo ±π (an interpretation in terms of a possible zero of the form factor was given in [11]). In this context, any alternative investigation of the scalar form factor based on analyticity, which might improve the knowledge of the scalar radius, is of great interest.
In a recent paper [12], a precise determination of the charge radius of the pion was obtained in the frame of a mixed dispersion formalism, using as input the phase of the electromagnetic form factor in the elastic region and the modulus above the first inelastic threshold. Unfortunately, in the case of the scalar form factor (1) no data on the modulus above the KK threshold are available. One can obtain however an integral constraint on the modulus squared of Γ π (t) along the cut using a formalism proposed a long time ago by Meiman [13] and Okubo [14], which exploits the dispersion relations for suitable QCD correlators, combined with unitarity and the positivity of the spectral functions. This formalism has been applied in the context of QCD for the first time in [15], and afterwards in many papers [16]- [33], being in particular a valuable tool for obtaining model-independent constraints on weak semileptonic form factors.
The Meiman-Okubo formalism has been applied also to the light-quark scalar correlator in [34,35], where it was used as a mean to derive a lower bound on the sum of the light quark masses. In the present paper we revisit the analysis reported in [34,35] bringing several improvements and updates. Thus, while in these works Watson theorem (5) was implemented only below 0.5 GeV, now it can be imposed up to the first relevant inelastic threshold, set by the KK channel, taking advantage of the recent progress in the determination of the pion-pion phase shifts. We include also higher terms in the expansion (3), which will lead to a more general constraint on the sum of the light quark masses, the derivatives of the scalar form factor at the origin and the phase shift δ 0 0 . Finally, we use the most recent calculation of the scalar correlator in perturbative QCD, available now to O(α 4 s ) [36]. The motivation of revisiting this analysis is the fact that precise results for both the light-quark masses and the pion scalar radius are now available from lattice calculations (for a recent review and earlier references see [37]). An updated, more precise independent constraint on these quantities is therefore of interest.
In the next section we describe the mathematical formalism, in Sec. III we discuss the input used in the calculations and in Sec. IV we present our results. Sec. V contains a discussion of the results in comparison with previous determinations and our conclusions.

II. DERIVATION OF THE BOUNDS
We consider the scalar correlator [34] Ψ(q 2 ) = i dxe iq·x 0|T (S(x)S † (0))|0 , written in terms of the operator S(x) defined in (2). The function Ψ(q 2 ) satisfies a dispersion relation which requires two subtractions. Therefore, as in [34] we shall consider the second derivative of Ψ, which is expressed on the Euclidian axis, i.e., for Q 2 = −q 2 > 0, as At large Q 2 , Ψ (Q 2 ) is given by the QCD perturbative expansion in powers of the renormalized strong coupling α s with negligible power corrections [36]: where the quark masses and the strong coupling are evaluated at a fixed scale µ 2 .
On the other hand, on the timelike axis the spectral function ImΨ(t) can be expressed, using unitarity, in terms of the contributions of hadronic states. Keeping the lowest two-pion contributions and using the positivity of the spectral function one obtains [34] where Γ π (t) is the pion scalar form factor defined in (1). From (7) and (9) one obtains the inequality where we denoted t p = 4m 2 π . By applying standard techniques of complex analysis, the right hand side of (10) can be further bounded from below by definite expressions involving values of the form factor at points inside the analyticity domain or the coefficients of the Taylor expansion (3) at t = 0.
In order to derive the optimal lower bound on the rhs of (10) with the constraints (3) and (5), we first perform the conformal mapping which maps the complex t plane cut for t ≥ t p onto the unit disk |z| < 1 such thatz(0) = 0,z(t p ) = 1 and the upper (lower) edge of cut becomes the unit semicircle ζ = e iθ with θ > 0 ( θ < 0). Then the inequality (10) can be written in the equivalent form wheret(z) = 4t p z/(1 + z) 2 is the inverse of (11) and φ(z) is an analytic function in |z| < 1 defined as with β Q = 1 + Q 2 /t p . By construction, φ(z) in an "outer function" [38], i.e., it has no zeros in |z| < 1, and its modulus on the boundary |z| = 1 is given by We emphasize that the absence of zeros in φ(z) guarantees that the bounds derived below are optimal. We further define a new function g(z), analytic in |z| < 1, by Then (12) implies the following inequality, valid for any K ≥ 1, where the real coefficients g k are defined by the Taylor expansion From (15) and (17) it follows that each g k is a linear combination of the derivatives of Γ π (t) at t = 0 up to the order k.
We now improve (16) by taking into account the additional constraint (5). We apply a technique of functional optimization based on functional Lagrange multipliers, proposed for the first time in [39], generalized and applied further in many papers [24,25,29,33]. We write below the solution of this optimization problem in our case (for a proof see [24,29]). Let be the image on the unit circle in the z plane of the point 4m 2 K + i situated on the upper edge of the cut [the point 4m 2 K −i being mapped onto exp(−iθ in )]. Then one obtains the stronger condition [24,29] where the function λ(θ) is the solution of the integral equation valid for θ ∈ (−θ in , θ in ), where the kernel is defined as in terms of the function We mention that integral equations of the form (20) are often encountered in solving functional optimization problems with constraints on the boundary implemented by Lagrange multipliers. The function λ (which is actually a generalized Lagrange multiplier) is a smooth, odd function of θ defined on the range (−θ in , θ in ), with values depending on the free parameters (r π , c π ) used in the optimization procedure. If the function Φ(θ) defined in (22) is sufficiently smooth, the integral equation (20) is of Fredholm type and is solved numerically by approximating it by a linear system of equations obtained by discretizing the integral. In the calculations performed in this work, the results proved to be very stable when the number of integration points was increased up to several hundreds.
The inequality (19), with Ψ (Q 2 ), given by (8), provides a model-independent relation between the sum of the light quark masses, the coefficients of the Taylor expansion of the scalar form factor, and its phase on the elastic part of the unitarity cut. In Sec. IV we shall explore numerically the consequences of this result. Before this we review in the next section the quantities used as input in the calculations.

III. INPUT
The perturbation expansion (8) of the scalar correlator is known at present in the MS scheme to O(α 4 s ) [36]. We used the coefficientsd 0,j written in [36,40] for 1 ≤ j ≤ 4 as polynomials of degree j of the quantity L Q = log Q 2 /µ 2 .  [48] and [49].
We have to specify the spacelike value Q 2 > 0 to be taken in the dispersion relation (7). The obvious require-ment is that the perturbative QCD expansion (8) must be meaningful. As discussed in [34] and in other applications of this formalism [24,26,29,30], the choice Q = 2 GeV is very reasonable for correlators involving light-quark operators, so we make this choice here. For illustration, as in [34], we take also Q = 1.5 GeV, which is still large enough compared to Λ QCD ∼ 0.300 GeV. Other choices of Q will be discussed in the last section. For the renormalization scale we made the choices µ = 2 GeV and µ = Q.
We obtained the strong coupling α s (µ 2 ) by using as input α s (m τ ) = 0.330 ± 0.010 [41] and evolving it to the scale µ by the renormalization group equation with β function to the same accuracy as the correlator [42,43]. The quark masses at µ = Q have been calculated starting from µ = 2 GeV and evolving them with the running to four loop from [44]. The apparent convergence of the expansion of Ψ (Q 2 ) is quite good: for Q = 2 GeV the contribution of the higher corrections in (8) is of 10% for α 2 s terms, 5% for α 3 s and 3% for α 4 s . For Q = 1.5 GeV the corrections are larger by a factor of about two 1 .
Considerable progress in the calculation of the ππ phase shifts using experimental data and Roy equations [45] has been achieved recently. In particular, several determinations of δ 0 0 (t) have been performed [46][47][48][49]. We used in the calculations the two solutions given in Eq. (11) and the Appendix of [48] and the most precise solution, with CFD parameters, specified in Eq. (A.3) and Table V of [49]. As can be seen from Fig. 1, where the central values of these phases are shown for √ t < 2m K , there are some differences between them, the prediction from [49] exhibiting slightly larger values around 0.8 GeV and a more pronounced increase near the KK threshold. We mention that the phase-shift δ 0 0 determined in [47], with the boundary values given in Eq. (72) of that work, is quite close to the second solution of [48], both exhibiting in particular a more moderate increase near the opening of the KK channel.

IV. RESULTS
By first setting K = 2 in the inequalities (16) and (19) and fixing the normalization Γ π (0), we obtain from these inequalities a lower bound on the sum m u +m d in terms of the scalar radius. Alternatively, the inequalities impose constraints on r 2 π s for input values of quark masses. Since the right hand sides of (16) and (19) are quadratic convex functions of r 2 π s , they lead to allowed ranges for this quantity, situated between a lower and an upper bound.
In Fig. 2 we show the parts of the boundaries of the  (16) and (19) leading to lower bounds on mu + m d for fixed r 2 π s or to lower bounds on r 2 π s for fixed mu + m d . Solid red (blue) lines: curves obtained with input phase equal to Solution 1 from [48], for Q = 2 GeV (1.5 GeV). Dashed red (blue) lines: curves obtained without phase input, for Q = 2 GeV (1.5 GeV).
domains (16) and (19) leading to lower bounds on m u + m d for fixed r 2 π s , or lower bounds on r 2 π s for fixed m u + m d , for a range of r 2 π s of physical interest. The renormalization scale has been fixed at µ = 2 GeV. The allowed values of the variables are situated above the curves. The parts of the boundaries giving upper bounds on r 2 π s involve values of this quantity too large to be of interest and are not shown in this figure (they will appear however in Figs. 3, 4 shown below).
The comparison of the dashed and solid curves, obtained without and with phase input, respectively, prove the significant improvement brought by the implementation of Watson theorem. Fig 2 shows also that the bounds obtained with Q = 1.5 GeV (blue curves) are more stringent than those obtained with Q = 2 GeV (red curves).
The results shown in Fig. 2 have been obtained with the central value Γ π (0) = 0.99 m 2 π and the first solution for the phase-shift δ 0 0 given in [48]. The phase shift from [49] leads to very close and slightly weaker results, while the second solution given in [48] leads to somewhat better (higher) bounds. It appears that a phase-shift δ 0 0 with a moderate increase near the KK threshold, such as exhibited by solution 2 of [48] and the phase-shift calculated in [47], leads to stronger lower bounds. These constraints can be used for testing the consistency of specific values for the sum of light quark masses and the pion scalar radius.
The study of the light quarks masses has a long history in the frame of low-energy effective theory of QCD [50,51]. The quantity difficult to estimate is the difference m u − m d , or the ratio m u /m d , for which recently an accurate value was obtained through the dispersive analysis [52] of the isospin-breaking decay η → 3π. We concentrate in our discussion on the recent lattice calculations of the quark masses and r 2 π s , summarized in the review [37]. From Table 7 of [37] one can obtain 13 lattice predictions for the sum m u + m d at the scale µ = 2 GeV, ranging from (5.60 ± 0.55) MeV [53,54] to (8.51±0.51) MeV [55], while four values of r 2 π s are given in Table 22 of the same review.
Using as input the lowest of these values, r 2 π s = (0.481 ± 0.062) fm 2 , obtained by the HPQCD Collaboration [56], and varying the normalization in the range (4) we obtained with the phase 2 from [48] a lower bound on m u +m d in the range (5.99−6.86) MeV for Q = 2 GeV, and in the range (6.76−7.74) MeV for Q = 1.5 GeV. The range obtained with Q = 1.5 GeV is in slight tension with the lowest value of m u + m d in Table 7 of [37], quoted above, obtained by MILC Collaboration [53,54].
For the other values of r 2 π s listed in Table 22 of [37], which are larger than the HPQCD prediction, one obtains from Fig. 2 smaller lower bounds on m u + m d , which are in no conflict with the values given in Table  7 of [37]. Adopting the central value r 2 π s = 0.61 fm 2 of the lattice calculations in [57][58][59] and using as input for the phase the solution 2 from [48], the most conservative lower bound, obtained with Q = 2 GeV and Γ π (0) = (0.99 ± 0.02) m 2 π , is This bound is consistent with all the lattice values listed in Table 7 of [37] and the average m u + m d = (6.85 ± 0.25) MeV quoted in the Review of Particle Physics [41]. We can generalize the constraints by setting in (16) and (19) the parameter K = 3, i.e., including also a curvature term in the expansion (3). Then, for a fixed value of the sum of quark masses, the inequalities (16) and (19) describe ellipses in the plane r 2 π s − c π s . We show for illustration in Fig. 3 the ellipses obtained using as input for the phase the solution 1 from [48] and the average value m u + m d = 6.85 MeV, quoted in [41]. The CFD parametrization [49] of the phase-shift leads to slightly larger domains, while for the solution 2 of δ 0 0 given in [48] the allowed domains are slightly smaller.
The comparison of the large ellipses with the small ones shows again the considerable effect of the implementation of Watson theorem. As before, the choice Q = 1.5 GeV leads to stronger constraints. In Fig. 4 we show for comparison the small ellipses, obtained with implementation of Watson theorem, for Q = 1.5 GeV and Q = 2 GeV.
From these ellipses one can read the values of the upper and lower bounds on r 2 π s obtained by using as input the sum of the quark masses. The upper bounds turn out to be very weak, so they are not useful for improving the accuracy of the dispersive predictions. On the other hand, as it was clear already from the previous discussion, the lower bounds are nontrivial. For instance, by varying the normalization (4) we obtain the conservative lower bounds r 2 π s ≥ 0.37 fm 2 for Q = 2 GeV and r 2 π s ≥ 0.53 fm 2 for Q = 1.5 GeV.
Another interesting property illustrated in Figs. 3, 4 is the strong correlation between the radius and the curvature at t = 0. In particular, for r 2 π s = 0.61 fm 2 , we obtained for the curvature the allowed range where we included the small variations due to the uncertainties of the phase δ 0 0 below the KK threshold, the normalization (4), the strong coupling α s and the truncation of the QCD expansion (8).
The allowed values given in (24) are higher than the prediction c π s ≈ 10 GeV −4 of the coupled-channel dispersive formalism [2,7]. A large value of the curvature would lead to rather large values of the relevant low energy constants (LEC) of χPT both for the two and three-flavour case. The study of these implications is beyond the scope of the present paper. It is important however to identify possible sources of systematic uncertainties that can affect both the dispersive approach and the bounds derived in this paper. We shall briefly discuss this problem in the next section.

V. DISCUSSION AND CONCLUSIONS
In this work, we revisited the application of the Meiman-Okubo formalism to the light-quark scalar correlator, improving the previous analysis [34] by the implementation of Watson theorem (5) up to the KK threshold and by the inclusion of higher derivatives in the Taylor expansion (3) of the pion scalar form factor. Moreover, recent progress in the determination of the scalar correlator in perturbative QCD and the phase shift δ 0 0 of pion-pion scattering has been taken into account.
Our result (19) is a model-independent constraint relating the sum of the light quark masses m u +m d (appearing in the expression (8) of the lhs) to the derivatives of the pion scalar form factor (entering the real coefficients g k ), and the phase shift δ 0 0 (t) below the KK inelastic threshold (appearing in the expression (22 of Φ(θ)). The phase of the form factor above the elastic region is not required in this approach. The results shown in Figs. 1 -3 illustrate the significant improvement brought by the implementation of Watson theorem along the entire elastic region up to the KK threshold 2 .
We have applied the above constraint for testing the consistency of the recent lattice results on the quark masses and the pion scalar radius r 2 π s . We found that, except a slight tension between the lowest value of r 2 π s given in Table 22 of the review [37] and the lowest value of m u + m d in Table 7 of the same review, the recent lattice determinations of the light quark masses and the pion scalar radius satisfy the consistency test.
As illustrated in Fig. 4, the sum of the light-quark masses and the phase below the inelastic threshold impose nontrivial constraints on the higher coefficients of the Taylor expansion (3) of the pion scalar form factor. The upper bounds on the scalar radius are very weak, being not useful for increasing the precision of the dispersive calculations of r 2 π s . On the other hand, the lower bounds turn out to be at the edge of the currently accepted values. Figure 4 shows also a strong correlation between the radius r 2 π s and the curvature c π s . This result is not sur-prising: strong correlations between the higher derivatives have been obtained also for other form factors in the frame of Meiman-Okubo formalism [23,25,33]. In the present case, somewhat surprising is the fact that the allowed range (24) for the curvature corresponding to r 2 π s = 0.61 fm 2 is considerably higher than the predictions of the dispersive treatment [2,7]. It is of interest therefore to discuss the possible systematic uncertainties that can affect the dispersive approach and the bounds derived in this paper.
The dispersive approach exploits unitarity, which relates the form factors to the meson-meson scattering amplitudes through a set of coupled homogeneous integral equations [2,6,7]. For solving these integral equations, each form factor is parametrized most generally as a polynomial P (t) multiplied by an Omnès function Ω(t), defined in terms of the phase δ(t) on the cut by In particular, for the form factor Γ π (t) the polynomial was taken as a constant, which implies that the form factor was assumed to have no zeros 3 . This assumption may be too restrictive: the presence of a polynomial which multiplies the Omnès function can lead to different predictions outside the limited interval of the unitarity cut where the coupled-channel equations are solved.
Another source of systematic uncertainty is the fact that the phase δ(t) is not known at higher energies, and some model-dependent assumptions about its behavior are required in the dispersive approach. One can check that the modulus of the Omnès function (25) behaves as t − δ(∞) π at large t. Therefore, if Ω(t) is multiplied by a polynomial, the asymptotic phase δ(∞) should be increased in order to ensure the asymptotic decrease as 1/t of the form factor, predicted by perturbative QCD. Although the higher derivatives at t = 0 are less sensitive to the phase at high energies, an anomalously large contribution of the inelastic channels and the presence of one or more zeros in the complex plane may have a significant contribution to the curvature.
On the other hand, the bounds derived in the present paper have been obtained with no assumptions about the phase above the KK threshold or the analytic expression of the form factor (i.e. the presence or absence of zeros). Below (24) we mentioned the small uncertainties due to the various pieces of the input. We shall now consider in more detail the normalization condition (4) and the perturbative QCD input as possible sources of systematic uncertainties.
As already mentioned, the error given in (4) might be underestimated. It is of interest to see how much would this uncertainty have to change in order to have a substantial impact on the results. We found that by increasing the error quoted in (4) by a factor of 10, the lower bound (24) on the curvature decreases slightly, becoming 31 GeV −4 . It turns out that in order to decrease the lower bound on the curvature below 20 GeV −4 , a huge increase of the error by a factor of about 40 would be required. We emphasize that the value of the form factor at the origin does not appear in the dispersive treatment [2,7], which involves only the ratio Γ π (t)/Γ π (0).
Turning to the perturbative QCD input, we note that from the maximum modulus principle it follows that the bounds depend in a monotonous way on the value of the lhs of the inequality (10), in the sense that larger values of Ψ (Q 2 ) for a fixed Q lead to weaker bounds. Making the very conservative assumption that the higher perturbative terms increase by a factor of 2 the value calculated from the sum in (8) truncated after four terms, we obtained the slightly larger interval c π s ∈ (31−36.7) GeV −4 . The bounds become weaker also if the value of the spacelike energy Q is increased (this was already illustrated in Fig. 4). Therefore, in order to obtain bounds of interest one should take the lowest Q for which the perturbative expansion is considered to be reliable. The range (24) was obtained by assuming that this value is Q = 2 GeV, which, as discussed in Sec. III, is reasonable for light-quark correlators evaluated on the spacelike axis. Assuming this value to be Q = 4 GeV gives the larger interval c π s ∈ (29−36) GeV −4 , while the choice Q = 10 GeV leads to the even larger range c π s ∈ (17 − 46) GeV −4 . In these calculations we used both scales µ = 2 GeV and µ = Q, which lead to very close results, and took the most conservative bounds.
Our analysis indicates that the model-independent bounds can accommodate the low values predicted by the older dispersive calculations only if the corrections to the normalization of the form factor at t = 0 are much larger than quoted in (4), or the perturbative QCD regime for the scalar correlator is assumed to start at quite large energies. Whether these strong assumptions are necessary or not remains to be established.
The results derived in this paper can be improved in principle by including in the unitarity sum for the spectral function ImΨ(t), besides the ππ states, also the contribution of the KK states, nonzero for t ≥ 4m 2 K , which can be expressed in terms of the kaon scalar form factor Γ K (t). One can thus include information about Γ K (0) available in χPT, however, as discussed in [6], the knowledge of this quantity is not very precise. In addition, one has to evaluate also the contribution of the unphysical cut of Γ K (t) along (4m 2 π , 4m 2 K ), which requires some model-dependent assumptions. Therefore, we shall not pursue this line here 4 .
We finally note that the uncertainties of the input quantities can be accounted for in a more realistic way by merging the present formalism with Monte Carlo simulations, as done in the recent analysis [12] of the pion charge radius.