A comprehensive revisit of the $\rho$ meson with improved Monte-Carlo based QCD sum rules

We improve the Monte-Carlo based QCD sum rules by introducing the rigorous H\"older-inequality-determined sum rule window and a Breit-Wigner type parametrization for the phenomenological spectral function. In this improved sum rule analysis methodology, the sum rule analysis window can be determined without any assumptions on OPE convergence or the QCD continuum.Therefore an unbiased prediction can be obtained for the phenomenological parameters (the hadronic mass and width etc.). We test the new approach in the $\rho$ meson channel with re-examination and inclusion of $\alpha_s$ corrections to dimension-4 condensates in the OPE. We obtain results highly consistent with experimental values. We also discuss the possible extension of this method to some other channels.


I. INTRODUCTION
QCD sum rule (QCDSR) is an important nonperturbative method in hadronic physics, introduced by Shifman et al. (SVZ) [1,2]. Due to its QCD based nature and minimally model dependent property, this semi-analytic approach has become a powerful weapon in the extensive studies on phenomenological hadronic properties. 1 Although fruitful results have been achieved within the framework of QCD sum rules, doubts on the predictive capability of this approach have never disappeared since in the original Laplace/Borel sum rules, the two external free parameters, i.e., Borel parameters τ and continuum thresholds s 0 , cannot be rigorously constrained. Intuitive choices of ranges for (τ , s 0 ) cannot completely exclude subjective factors in determination of the phenomenological outputs.
Two main directions have been developed to address the above issues. First, the sum rule stability criteria [4,5] in Laplace sum rules (LSR), systematically developed by Narison, have been widely and successfully used to study QCD phenomenological properties. 2 In the widely used "single narrow resonance minimal duality ansatz" where the hadronic widths are not involved, the mass curves present stability (if any) in τ or s 0 at the extremums, where the mass depends weakly on (τ , s 0 ) and thus optimal predictions can be achieved if validity of operator product expansion (OPE) truncation and pole contribution dominance are simultaneously ensured. In the cases where the s 0 -stability is absent, the finite energy sum rules (FESR) [12] can complement the systematics of the method [11,13,14]. These approaches require an assessment of whether the optimized values for τ and s 0 are in the sum rule region of validity.
The Monte-Carlo based weighted-least-squares fitting procedure is another quite successful approach to allow reliable sum rule predictions [15][16][17][18][19][20][21][22]. As originally suggested by SVZ [1,2], the sum rules should be analyzed in a certain range of the Borel parameter, called the "sum rule window", in which the OPE series reach convergence and the continuum contribution are suppressed. Directly following this discussion, Leinweber introduced the Monte-Carlo based weightedleast-squares fitting procedure [15] in the sum rule numerical analysis. Within this approach, phenomenological outputs, such as the mass, the decay constant, and the continuum threshold etc., can be obtained by minimizing the weighted-χ 2 function between the OPE and phenomenological expressions of the correlation function within the sum rule window. This method has some obvious advantages: • the continuum threshold s 0 can be determined from the fitting procedure; • different parametrization models for the hadronic spectral function can be dealt with therefore it is possible to obtain predictions both for hadronic mass and width with an appropriate spectral density; • uncertainty analyses and dependence of the outputs on the phenomenological input parameters can be provided via the Monte-Carlo based numerical procedure.
However, as in the original (as well as many common uses of) QCD sum rules, the sum rule window cannot be rigorously constrained. 3 The usual imposed conditions, the pole contribution > 50% of total phenomenological expression and the highest dimension operator (HDO) contribution < 10% of the total OPE are based upon assumptions, and cannot always be ensured especially for multi-quark state sum rules where the OPE convergence is slower, which makes artificial adjustments of these constraints inevitable.
The Hölder inequalities for QCD sum rules, obtained from the requirement that the imaginary part of a correlation function should be positive because of its relation to a hadronic spectral function (e.g., physical cross sections for electromagnetic currents), were first introduced as fundamental constraints on the sum rules by Benmerrouche et al. [23], then were used in many research works on QCD sum rules, including constraints on the QCD sum rule window [24][25][26][27][28]. However, these inequalities have mainly been applied in analyses of the simplified "single narrow resonance minimal duality ansatz".
In fact, most hadrons listed in Review of Particle Physics [29] are resonances rather than stable particles, thus the mass and width are equally important in describing the properties of these hadrons, and equal attention should be paid to them. Some researchers replace the delta function in spectral density directly with a Breit-Wigner (BW) type function in order to obtain information for the hadron width in QCD sum rules [19,30,31]. In this paper, we will introduce a BW type parametrization for the phenomenological spectral function that satisfies the low-energy theorem for the form factor, and apply the rigorous Hölder-inequality-determined sum rule window in the Monte-Carlo based fitting procedure. We will reanalyze the phenomenological properties of the ρ meson with this new systematic numerical procedure. A comprehensive study will be presented, including the explicit reexamination of some α s corrections in the OPE 4 and uncertainty analysis for the fitted results. Extensions of our analysis to include a smooth transition to the QCD continuum are also presented. We will conclude the paper by summarizing our calculation and analysis and discussing the possible extension of our approach.

II. OPERATOR PRODUCT EXPANSION FOR CORRELATION FUNCTION FOR VECTOR I = 1 CURRENT
The starting point of QCD sum rules is to calculate the correlation function for a specific current. In this paper, we consider the light vector current with isospin 1, i.e., j µ =ūγ µ d (and 1 √ 2 (ūγ µ u −dγ µ d),dγ µ u), which has the same quantum numbers as the ρ meson. Because of conservation of the vector current, the correlation function for current j µ can be written as where Π(q 2 ) is a Lorentz scalar function which can be calculated theoretically using operator product expansion methods [33]. The well-confirmed and widely used QCD expression of Π(q 2 ), including leading-order (LO) and nextto-leading order (NLO) perturbative terms and LO non-perturbative terms up to dimension-6 condensates reads [1,2]: where κ is the factorization violation factor which parameterizes the deviation of the four-quark condensate from a product of two-quark condensates. In addition to these terms, some researchers have also calculated the NLO corrections in α s to dimension-4 condensates [32,34], which however are not included in the existing sum rule numerical analyses for ρ meson [15,30]. Although the results in [32,34] suggest these corrections are relatively small, they can play an important role in determining the sum rule window from the Hölder inequalities, motivating our independent calculation with a different method. Thus, we calculate these contributions using the external field method in Feynman gauge 5 [33,37] to 3 A remedial strategy is to determine to optimal outputs by observing the variation of the outputs on the change of the conditions imposed on the sum rule window. If the outputs are not sensitive to the variation of the range of the window, reliable predictions can be obtained by demanding the OPE converge in a proper trend, as did in [22] for the exotic hybrids. 4 These radiative corrections to dimension-4 quark and gluon condensates are complicated calculations but may play an important role in determining the range of τ where the Hölder inequalities are satisfied. The previous results of these radiative corrections obtained using a projector method in [32] were not used in the existing numerical analyses [15,30], motivating our independent calculation. 5 Similar calculations can be seen in [35,36] for the (1 −+ ,0 ++ ) light hybrid states. examine the previous results obtained in [32,34]. Our results read which is consistent with the previous results in Ref. [32,34]. The whole calculation is much more complicated than the LO calculation, thus for clarity, we give all the related Feynman diagrams and explicit results of our calculation in the Appendix. In the following section, we will use to reanalyze the sum rules for the ρ meson, and discuss the effects of these α s corrections.

III. MONTE-CARLO BASED QCD SUM RULES FOR ρ MESON
To obtain a QCD sum rule which can be used to predict hadronic properties, we first need to Borel-transform the theoretical representation of the correlation function, i.e., Π (OPE) (q 2 ). After some calculations, we obtain whereB is the Borel transformation operator and α s (1/τ ) = 4π/(9 ln(1/(τ Λ 2 QCD ))) is the running coupling constant for three flavors at scale 1/ √ τ . We have considered the renormalization-group (RG) improvement of the sum rules [38] and anomalous dimensions for condensates [1,2] in Eq. (5), where µ 0 is the renormalization scale for condensates.
Eq. (5) provides the theoretical representation of Borel-transformed correlation function, but in order to obtain a QCD sum rule, we still need the phenomenological representation obtained by constructing a phenomenological spectral density model.
If we insert a complete set of one-particle states other states" into correlation function, we will reach a phenomenological spectral density with a delta function δ(s − m 2 ρ ), which is widely used in traditional QCD sum rules. However, the ground state of light vector I = 1 meson, i.e., ρ meson, is a resonance which is far away from a stable particle, thus it's more appropriate to insert two-particles intermediate states into the correlation function since the two-pion decay mode is the dominant decay mode for ρ meson. By inserting two-pion intermediate states into Eq. (1), e.g., inserting for correlation function of current j µ =dγ µ u, and using the Cutkosky's cutting rules [39], the phenomenological expression for ImΠ(s) can be found [40]: |F π (s)| 2 + contributions from excited states and continuum, where m π is the mass of pion, and has been used, where F π (s) is the electromagnetic form factor which is normalized as F π (0) = 1. Furthermore, since the main contribution to F π (s) comes from the ρ meson, F π (s) should have a pole at s = m 2 ρ − im ρ Γ ρ , where m ρ and Γ ρ are the mass and width of ρ meson respectively. Following Ref. [41] we use a Breit-Wigner form function to construct a model for the form factor as follows which meets the above requirements (including the low-energy theorem F π (0) = 1), and is more simple than the Gounaris-Sakurai parameterized form factor used in Ref. [40] and the parameterized form factor used in Ref. [42,43].
As outlined below, the phenomenological spectral function normalization constrained by a low-energy theorem enters our analysis in a meaningful way because we work directly with the Borel-transformed correlation function R(τ ) rather than the typical sum rule approach which uses normalization-independent ratios of sum rules. For the excited states and continuum (ESC) contributions in the spectral density, we still use the same model as traditional QCD sum rules, i.e., we use a spectral density as follows in this paper: where s 0 is the continuum threshold separating the contributions from excited states, and we have omitted the small mass of the pion. By using the dispersion relation, we can obtain the phenomenological representation for Borel-transformed correlation function, which has a following form: Then the master equation for QCD sum rule can be obtained by demanding the equivalence between Eq. (5) and (10): which can be used to obtain the predictions for s 0 , m ρ and Γ ρ . Obviously, because of the truncation of OPE and the simplicity of the phenomenological spectral density, Eq. (11) can not be valid for all τ , thus one requires a sum rule window in which the validity of the master equation can be established. Usually, a sum rule window is determined by demanding the highest dimension operator contributions be limited to 10% of total OPE contributions and the continuum contributions be limited to 50% of total phenomenological contributions [15,19,30]. However, the setting of percentage 10% and 50% is somewhat arbitrary.
Benmerrouche et al. presented a method based on the Hölder inequality which provides fundamental constraints on QCD sum rules [23]. The Hölder inequality for integrals defined over a measure dµ [44,45] is where 1/p + 1/q = 1 and p, q ≥ 1. Since ImΠ(s) is positive because of its relation to spectral functions (and in some cases to physical cross sections), it can serve as the measure dµ = ImΠ(s)ds in Eq. (12). By placing the excited states and continuum contributions on the OPE side, we obtain then the Hölder inequality for QCD sum rules can be written as [23] where f (s) = e −ωsτ1 and g(s) = e −(1−ω)sτ2 were used. The parameter ω is defined as ω = 1/p, which satisfies 0 ≤ ω ≤ 1. For parameters τ 1 and τ 2 we demand τ 1 < τ 2 . Following Ref. [23], we will perform a local analysis on Eq. (14) with τ 2 − τ 1 = δτ = 0.01 GeV −2 . The only starting point of this inequality is that ImΠ(s) should be positive because of its relation to physical spectral functions (including cross sections in the case of electromagnetic currents), thus Eq. (14) must be satisfied if sum rules are to consistently describe integrated physical spectral functions. The Hölder inequality can be used to check if a sum rule window is reliable [23,24], or to give extra constraints on some physical quantities [25]. It can also be used to determinate the QCD sum rule window [26][27][28]. In this paper, we will use Eq. (14) as the only constraint on the QCD sum rule window: by choosing the maximally allowed region of the Hölder inequality, we will determinate a sum rule window directly. In order to match the two sides of master equation (11) in the sum rule window, a weighted-least-squares method [15] will be used in this paper. By randomly generating 200 sets of Gaussian distributed phenomenological input parameters with given uncertainties at τ j = τ min + (τ max − τ min ) × (j − 1)/(n B − 1), where n B = 21, we can estimate the standard deviation σ OPE (τ j ) for R (OPE) (τ j ). Then, the phenomenological output parameters s 0 , m ρ and Γ ρ can be obtained by minimizing After obtaining σ OPE , 2000 sets of Gaussian distributed input parameters with same given uncertainties will be generated, we will minimize χ 2 to obtain a set of phenomenological output parameters for each set and the uncertainties of s 0 , m ρ and Γ ρ can be estimated based on these results.

IV. NUMERICAL RESULTS
In our numerical analysis, we use the central values of input parameters (at µ 0 = 1 GeV) from a recent review article of QCD sum rules [5]. These values read where the factor κ indicates the violation of factorization hypothesis in estimating the dimension-6 quark condensates. The size of κ have been observed in different channels to be 2-4 [46][47][48], which is the scale we will consider in our analysis.
Before fitting the two sides of the master equation (11), we should determinate the sum rule window first. In FIG. 1(a) and 1(b), we plot the region which is allowed by the Hölder inequality (14) with κ = 2, 3 respectively. For larger factorization violation factor κ, the allowed region shrinks, however, the main characteristics of the allowed region, such as the existence of a lower bound on s 0 , and the bound on τ , will persist for any value of κ from 2 to 4. From FIG. 1(a) and 1(b), we find that the α s corrections to α s G 2 and m qq q extend the allowed region to a higher τ region and lower s 0 region. This extension is more obvious in the case with a small κ. Thus according to the Hölder inequality, the α s corrections to dimension-4 operator condensates extend the validity of the sum rule master equation. As a comparison, we also plot the region obtained by demanding HDO/OPE < 10% and ESC/total < 50% with s 0 > 1 GeV 2 . However, because of the small order of magnitude, the α s corrections to dimension-4 operators seem to not have important effects in the determination of the sum rule window from 50%-10% method. Whether these corrections are included or not, the allowed region does not change a lot, thus we do not plot the region for which such corrections are not considered. From FIG. 1, we also observe that the constraint on the upper bound of τ from the Hölder inequality is more strict than the 10% method for s 0 providing that s 0 is not too large. However, the constraint on the lower bound of τ is looser than that from the 50% method, the largest proportion of ESC can be 60%-65% of total phenomenological contributions.  (14) for κ = 2 (a) and κ = 3 (b). The region with (blue) dot or (red) line is allowed for sum rule with or without αs corrections to αsG 2 and mqqq respectively. The region surrounded by the (black) dashed line is obtained by the requirement of HDO/OPE < 10% and ESC/total < 50% with αs corrections to αsG 2 and mqqq .
We can not obtain a fixed sum rule window from FIG. 1 because s 0 is not determined at first, thus in practice, we will initially "guess" an s 0 , then we will obtain a sum rule window [τ min , τ max ] from FIG. 1, where τ min and τ max are respectively the lower bound and upper bound of the allowed τ region with the guessed s 0 . After finishing the analysis by using this sum rule window, we can check our initial choice on s 0 with its fitted value and adjust it iteratively until it is consistent with the results of the analysis. After finishing this iterative procedure, the sum rule window is completely determined, all outputs will be obtained self-consistently from the final sum rule window.
include Π α d=4 χ 2 /10 −4 s0/GeV 2 mρ/GeV Γρ/GeV  The uncertainties of all input parameters in (16) are set to 10% which is a typical uncertainty in QCDSR [15]. After finishing the fitting procedure, we obtain the results which are listed in TABLE I, where the median and the asymmetric standard deviations from the median for all output parameters are reported. All results are the final statistical results from 2000 fitting samples, and we report both the results with and without α s corrections to α s G 2 and m qq q .
From TABLE I we find that the uncertainties of phenomenological output parameters are less than 10% for κ = 3, 4. Only with κ = 2, the uncertainty of Γ ρ will reach to about 15%, however, even in this case, the uncertainties of s 0 and m ρ are still less than 10%. These results imply that the fitted results are very stable with different input parameters. After including the α s corrections, we find all values of output parameters reduce about 4%-6% for κ = 4, for small κ, the reduction is even more apparent. This result is in contrast with the result in Ref. [32] where the α s corrections increased the predicted value of the ρ meson mass.
We also find besides extending the allowed region by the Hölder inequality and reducing the output parameters, the α s correction also reduce the value of χ 2 . For a least-squares fitting, a smaller χ 2 means a better fitting, thus the fitted results are more reliable if these α s corrections to dimension-4 operators are included. Based on the 2000 fitting samples, the correlations of all parameters can also be obtained. From the scatter plot FIG. 2(a), we find that there exists very strong positive correlation between m ρ and κα s qq 2 , which means a larger contribution from dimension-6 operator condensate will lead to a larger m ρ . From FIG. 2(b), we also learn that a larger m ρ is accompanied with a larger Γ ρ , thus the determination of an exact value of dimension-6 operator condensate is extremely important. We also find that there exists weak negative correlation between α s G 2 and all output parameters, thus any correction to α s G 2 is important. However, because of the small order of magnitude, m qq q is not important in the sum rules for the ρ meson. Finally, a larger Λ QCD also leads to smaller output parameters, there exists weak negative correlation between them.
From TABLE I we find if we choose κ ∼ 3, the fitted result (with 10% uncertainties on input parameters) will cover the physical mass of ρ [29]. If we want to find an exact match between the predicted mass and the physical mass, then κ = 2.8 is the best choice, which leads to the result as follows s 0 = 1.71 +0. 12 −0.13 GeV 2 , m ρ = 0.774 +0.037 −0.046 GeV, Γ ρ = 0.126 +0.010 −0.012 GeV. (17) However, the predicted width is lower than the physical value [29] by 16%, i.e., we meet a similar problem as Ref. [40]. However, we can conclude from the scatter plot between m ρ and Γ ρ that it can not be explained by the uncertainty in the value of the four-quark condensate because if the width matches the experimental value, then the mass will be too large. Shuryak once parameterized the spectral density from the experimental data for the vector I = 1 channel by the following function [49] 1 π ImΠ ρ (s) = 3 2π 2 where E 0 = 1.3 GeV, δ = 0.2 GeV andᾱ s (s) = 0.7/ ln( √ s/0.2 GeV). This function includes a resonance peak (the first term of Eq. (18)) and excited states and continuum contributions (the second term of Eq. (18)) which smoothly transition to the perturbative result of the spectral density as s increases. Now we can compare our spectral density model with Eq. (18). From FIG. 3, we can conclude that the simple spectral density we are using has characterized the main behavior of ρ resonance peak. However, we also notice that the region around s 0 is the main incompatible region in the spectral density, although the value of s 0 with κ = 2.8 is very close to E 2 0 , the "step" at s 0 in the spectral density is unnatural. A smoothly varying spectral density should be used in this region, and we conjecture that a better description of the ESC may lead to a better prediction for Γ ρ . To investigate our conjecture, we replace the θ function in the traditional ESC model 1 π ImΠ (OPE) (s)θ(s − s 0 ) with a smooth function 6 to construct a new ESC model for the phenomenological spectral density as follows where we still associate s 0 with the continuum threshold while λ (≥ 0 GeV 2 ) as the parameter which controls the transition of 1 π ImΠ (ESC) (s) to the QCD continuum. Obviously, when λ → 0 GeV 2 , the new ESC model will recover the same ESC as in Eq. (9). Conversely, a larger λ will lead to a flatter ESC contribution in the phenomenological spectral density.
To obtain an impression of how the fit is affected by λ, we first fix the sum rule window from the λ = 0 Hölder inequality constraint, then perform the least-χ 2 fit (15) for selected values of λ. Because all fits (with same κ) use the same sum rule window, we can identify the preferred range of λ from values of χ 2 . In Fig. 4, we plot this minimized χ 2 as a function of λ for κ = 3.0, 3.5, 4.0. Although the evidence for a non-zero optimum value of λ is marginal, the rapid increase in χ 2 for large λ provides an upper bound of approximately λ 1 GeV 2 . By setting λ = 0.2, 0.4, 0.6, 0.8, 1.0 GeV 2 , we can obtain the allowed (τ , s 0 ) regions from the Hölder inequality respectively, then the iterative fitting procedure introduced in the previous section can be invoked. We plot ratios of the fitted mass and width to the experimental values as λ increases in Fig. 5 (related detailed fitted results are listed in Table II). From Fig. 5 and Table II, we find that a larger value of λ will generally lead to a smaller m ρ , and a larger Γ ρ (the fitted Γ ρ is not sensitive with λ with λ < 0.4 GeV 2 ), furthermore, increasing κ, increases both m ρ and Γ ρ . These tendencies demonstrate that with an appropriately smoothed ESC and 0.4 GeV 2 < λ < 1 GeV 2 , we may obtain results for Γ ρ and m ρ statistically consistent with their experimental values.
To obtain an intuitive impression on the behavior of the ESC contributions, we plot the two spectral density models used in this paper in Fig. 6, where we also plot the ESC from Shuryak's spectral density (18). From this figure, we find that our new ESC model describes the main behavior of Shuryak's ESC [49] better than the traditional ESC model.
Finally, as a comparison, we also try the fitting procedure with the traditional phenomenological narrow width κ = 3.0 κ = 3.5 κ = 4.0 s0/GeV 2 mρ/GeV Γρ/GeV s0/GeV 2 mρ/GeV Γρ/GeV s0/GeV 2 mρ/GeV Γρ/GeV where f ρ is the coupling constant of ρ meson. However, the outputs of s 0 with κ = 2-4 are too small, thus the sum rule window departs from the allowed region from the Hölder inequality. Even without violation of factorization, we still can not find a result in this procedure which is completely consistent with the Hölder inequality.

V. SUMMARY AND DISCUSSION
In this paper, we have introduced a systematic sum rule approach to give predictions for the hadronic mass and width. This approach is a synthesis of the Monte-Carlo based QCD sum rule and the Hölder-inequality-determined sum rule window, accompanied with a phenomenological spectral density of Breit-Wigner form that is normalized by the low-energy theorem. We apply this approach to the ρ meson to re-analyze its phenomenological properties and the corresponding theoretical uncertainties. We also calculate the α s corrections to dimension-4 condensates in the OPE of light vector I = 1 current correlation function by using the external field method in Feynman gauge. Our calculation confirm the previous results obtained using the projector method [32], which have not previously been considered in the existing sum rule analyses. Considering these higher order effects in the perturbation series, we conducted a comprehensive reanalysis with the new methodological approach. The findings of our numerical analysis are: • The Breit-Wigner type spectral density model, which provides a better description of physical spectral density than the traditional single narrow resonance spectral model, plays an important role in our procedure. First, the sum rule window can be completely determined by the Hölder inequality in this spectral density model, thus we avoid the often criticized 50%-10% assumptions used to constrain the sum rule window in previous Monte-Carlo based QCD sum rules analysis [15,19,21]. Meanwhile, the traditional single narrow resonance phenomenological spectral density is too over-simplified to find a fitted result consistent with the Hölder inequality. Second, because of the low-energy theorem normalization for the electromagnetic form factor in our phenomenological spectral density, the degrees of freedom during the fitting procedure reduces to three phenomenological output parameters, i.e., s 0 , m ρ and Γ ρ and improves the ability of QCD sum rules to predict the width of ρ meson.
The three-parameter model is a clear improvement on the unphysical Γ ρ = 0 result [21] that is found in the four-parameter spectral density model 1 π ImΠ(s) = f 2 • Although most works in the literature do not consider the effect of α s corrections to dimension-4 operators, we find they play an important role in the present QCD sum rules analysis for ρ meson by extending the s 0 -τ region allowed by the Hölder inequality and improving the goodness for fit between the theoretical side and the phenomenological side of the master equation for QCD sum rules. These radiative corrections also reduce the values of s 0 , m ρ and Γ ρ by 4%-16%, depending on the value of κ. Thus it is important to include these corrections in our QCDSR approach, although the order of magnitude of these corrections is small.
• With 10% uncertainties for input parameters, the optimal phenomenological results in our numerical analysis are: m ρ = 0.774 +0.037 −0.046 GeV, Γ ρ = 0.126 +0.010 −0.012 GeV and s 0 = 1.71 +0.12 −0.13 GeV 2 with κ = 2.8 from the traditional ESC model. All the outputs are consistent with the experimental values for the ρ meson. The factorization violation factor κ ∼ 3 is also the expected size estimated from many different channels. Extending our phenomenological model to include a smooth transition to the QCD continuum (similar to Ref. [49]) leads to generally better agreement with the experimental values.
The successful application of our approach in the ρ channel greatly motivates us to extend it to other channels. Since similar phenomenological spectral densities can be constructed by inserting a complete set of intermediate states for other hadrons which have dominant decay modes, such extensions are expected to be applicable. For example, by inserting two-pion intermediate states in the correlation function for I = 0 scalar current, we can connect the phenomenological spectral density with scalar pion form factor F s (s). Then the result of F s (0) = m 2 π in chiral perturbation theory [50] can be used to reduce the number of the phenomenological parameters to three. Furthermore, the Hölder inequality constraint to determine the sum rule window of validity is similarly adaptable to other channels. Subsequent works can examine the validity of our new method and help us further improve it.