Global determination of two-meson distribution amplitudes from three-body $B$ decays in the perturbative QCD approach

We perform a global analysis of three-body charmless hadronic decays $B\to VP_3\to P_1P_2P_3$ in the perturbative QCD (PQCD) approach, where $V$ denotes an intermediate vector resonance, and $P_i$, $i=1,2,3$, denote final-state pseudoscalar mesons. Fitting the PQCD factorization formulas at leading order in the strong coupling $\alpha_s$ to measured branching ratios and direct $CP$ asymmetries, we determine the Gegenbauer moments in the two-meson distribution amplitudes (DAs) for the meson pairs $P_1P_2=\pi\pi, K\pi, KK$. The fitted Gegenbauer moments are then employed to make predictions for those observables, whose data are excluded in the fit due to larger experimental uncertainties. A general consistency between our predictions and data is achieved, which hints the validity of the PQCD formalism for the above three-body $B$ meson decays and the universality of the nonperturbative two-meson DAs. The obtained two-meson DAs can be applied to PQCD studies of other multi-body $B$ meson decays involving the same meson pairs. We also attempt to determine the dependence of the Gegenbauer moments on the meson-pair invariant mass, and demonstrate that this determination is promising, when data become more precise.

evolution equation [24][25][26][27], i.e., the Gegenbauer polynomials C 3/2 n (2x − 1), based on the conformal symmetry. This expansion follows that for a hadron DA exactly. As to the expansion in ζ, one employs the partial waves for the produced meson pair, i.e, the Legendre polynomials P l (2ζ − 1), noticing the relation 2ζ − 1 = cos θ with θ being the polar angle of a meson in the centerof-mass frame of the meson pair [28]. The expansion of a two-meson DA in terms of the two sets of orthogonal polynomials then reads [22] Φ(x, ζ, ω 2 ) = where B nl (ω 2 ) are the ω 2 -dependent coefficients, N c = 3 is the number of colors, and l = 0, 1, 2,... denote the S-wave , P -wave, D-wave,... components, respectively. The time-like form factor B 0l (ω 2 ), which normalizes each of the partial-wave component, contains both resonant and nonresonant contributions. Some form factors, such as the time-like pion form factor that receives contributions from the series of ρ resonances, have been constrained stringently by experimental data [29]. The other coefficients B nl (ω 2 ), referred to as the Gegenbauer moments, are still quite uncertain due to a lack of systematic nonperturbative studies. Note that these Gegenbauer moments differ from those in the DA for a specific resonance which strongly decays into the meson pair, because, as stated above, a two-meson DA collects contributions from a series of resonances as well as nonresonant contributions. Moreover, they are ω 2 -dependent, a feature dramatically distinct from the Gegenbauer moments for a meson DA. It has been observed [30] that the Gegenbauer moments of a P -wave di-pion DA differ from those of the ρ(770) meson DA. Therefore, it is essential to determine the Gegenbauer moments for two-meson DAs in order to improve the precision of theoretical predictions for multi-body B meson decays in factorization frameworks.
We will perform a global fit of the Gegenbauer moments in two-meson DAs to measured branching ratios and direct CP asymmetries in three-body charmless hadronic B meson decays B → V P 3 → P 1 P 2 P 3 in the PQCD approach, where V stands for an intermediate vector resonance, and P i , i = 1, 2, 3, stand for final-state pseudoscalar mesons. As the first attempt to a global determination of two-meson DAs, we focus on the P -wave components, and employ the LO PQCD factorization formulas for decay amplitudes. We establish a Gegenbauer-moment-independent database, by means of which each decay amplitude is expressed as a combination of the relevant Gegenbauer moments in the two-meson DAs. The Gegenbauer moments in the DAs for the mesons P 3 = π, K are input from the global analysis of two-body B meson decays in Ref. [31]. The leading-twist (twist-2) and next-to-leading-twist (twist-3) DAs for the pairs P 1 P 2 = ππ, Kπ and KK with the intermediate vector mesons V = ρ, K * and φ, respectively, are then fixed in the global fit. Because the current data for three-body B meson decays are not yet precise enough to determine the ω 2 dependence of the Gegenbauer moments, we first treat them as constant parameters defined at the initial scale 1 GeV. One or two Gegenbauer moments for each of the above two-meson DAs are obtained with satisfactory fit quality, depending on the abundance of available data. It is noticed that the results and the precision of the extracted two-meson DAs depend on the number of the Gegenbauer moments considered in the fit: when more Gegenbauer moments are introduced into the Kπ DAs, the quality of the fit is improved at the cost of amplified uncertainties for fit outcomes.
The determined Gegenbauer moments are then employed to make predictions for those observables, whose data are excluded in the fit due to larger experimental errors. A general consistency between our predictions and data for various modes is achieved, except those which suffer significant subleading corrections according to previous PQCD studies, such as the B 0 → π 0 (ρ 0 →)ππ decay [32,33]. The consistency hints the validity of the PQCD formalism for three-body hadronic B meson decays and the universality of the nonperturbative two-meson DAs. The ππ, Kπ and KK twist-2 and twist-3 DAs presented in this work are ready for applications to PQCD investigations of other multi-body B meson decays involving the same meson pairs. Our formalism can be extended to global fits for other two-meson DAs of various partial waves straightforwardly. It can be also generalized to include higher-order and/or higher-power corrections to PQCD factorization formulas [34], when they are available, so that more accurate two-meson DAs are attainable in a systematic way.
As a more ambitious attempt, we explore the dependence of the Gegenbauer moments in the di-pion DAs on the pion-pair invariant mass squared ω 2 . It is unlikely to determine the exact ω 2 dependence from current data, so we simply parametrize the Gegenbauer moments up to the first power in ω 2 , following their series expansion derived in Ref. [22]. The global fit shows that at least the linear term in one of the twist-3 di-pion DAs can be constrained effectively. It indicates that the determination of the ω 2 -dependent Gegenbauer moments in two-meson DAs is promising, when data become more precise in the future.
The rest of the paper is organized as follows. The kinematic variables for three-body hadronic B meson decays are defined in Sec. II, where the dependence on final-state meson masses is included to describe the phase space accurately. The parton kinematics and hard kernels are then refined, such that SU (3) symmetry breaking effects in the decays can be evaluated more precisely. The considered two-meson P -wave DAs are parametrized, whose normalization form factors are assumed to take the relativistic Breit-Wigner (RBW) model [35] or the Gounaris-Sakurai (GS) model [36]. We explain how to perform the global fit, present and discuss the numerical results, and try to extract the ω 2 dependence of the Gegenbauer moments in Sec. III, which is followed by the Conclusion. We collect the PQCD factorization formulas for the decay amplitudes with numerous refined hard kernels in the Appendix.

A. Kinematics
Consider the charmless B meson decay into three pseudoscalar mesons via a vector intermediate resonance, B(p B ) → V (p)P 3 (p 3 ) → P 1 (p 1 )P 2 (p 2 )P 3 (p 3 ), with the meson momenta p B = p + p 3 and p = p 1 + p 2 . We work in the B meson rest frame and parametrize the relevant momenta in the light-cone coordinates as where m B is the B meson mass, and k B , k and k 3 are the valence quark momenta in the B meson, the meson pair, and the bachelor meson with the parton momentum fraction (transverse momenta) x B , z and x 3 (k BT , k T and k 3T ), respectively. That is, we have chosen the frame, such that the meson pair and the bachelor meson move in the directions n = (1, 0, 0 T ) and v = (0, 1, 0 T ), respectively. Since the parton momentum k (k 3 ) is aligned with the meson pair (bachelor meson), its small minus (plus) component has been neglected. We have also dropped the plus component k + B , because it does not appear in the hard kernels for dominant factorizable contributions. In the above expressions, the functions f ± and g ± are written as with the ratio r 3 = m 2 P3 /m 2 B (s) and η = ω 2 /m 2 B (s) , m P3 being the bachelor meson mass and ω 2 = p 2 being the invariant mass squared of the meson pair. For a P -wave meson pair, we introduce the longitudinal polarization vector We derive the meson momenta p 1 and p 2 , from the relation p = p 1 + p 2 and the on-shell conditions p 2 i = m 2 Pi , i = 1, 2, with the mass ratios r 1,2 = m 2 P1,P2 /m 2 B . The variable ζ + (r 1 − r 2 )/(2η) = p + 1 /p + bears the meaning of the meson momentum fraction up to corrections from the final-state meson masses. Alternatively, one can define the polar angle θ of the meson P 1 in the P 1 P 2 pair rest frame. The transformation between the B meson rest frame and the meson pair rest frame leads to the relation between the meson momentum fraction ζ and the polar angle θ, with the bounds We emphasize that the parametrization with the exact dependence on the final-state meson masses in Eq. (5) is crucial for establishing Eq. (6), such that the Legendre polynomials in Eq. (1) correspond to the partial waves of the meson pair exactly.
The branching ratio for a three-body B meson decay is given by [37] with the B meson lifetime τ B . The direct CP asymmetry A CP is then defined as The decay amplitude A, according to the factorization theorem stated in the Introduction, is expressed as where Φ B (Φ P3 ) is the B (bachelor) meson DA, and the two-meson DA Φ P1P2 absorbs the nonperturbative dynamics in the production of the meson pair P 1 P 2 . The symbol ⊗ denotes the convolution of the above factors in parton momenta.

B. Distribution Amplitudes
The light-cone hadronic matrix element for a B meson is parametrized as [38][39][40][41][42] where q represents a u, d or s quark. The two wave functions φ B andφ B in the above decomposition, related to φ + B and φ − B defined in the literature [43] It has been shown that the contribution fromφ B is of next-to-leading power and numerically suppressed [39,40,44], compared to the leading-power contribution from φ B . Taking the PQCD evaluation of the B → π transition form factor F B→π 0 in Ref. [44] as an example, we find that theφ B contribution to F B→π 0 is about 20% of the φ B one. The higher-twist B meson DAs have been systematically investigated in the heavy quark effective theory [45], which are decomposed according to definite twist and conformal spin assignments up to twist 6. In principle, all the next-to-leading-power sources should be included for a consistent and complete analysis, which, however, goes beyond the scope of the present formalism. Therefore, we focus only on the leading-power component with the impact parameter b being conjugate to the parton transverse momentum k BT . The B meson DA φ B (x, b) is chosen as the model form widely adopted in the PQCD approach [38][39][40][41][42]46], where the constant N B is related to the B meson decay constant f B through the normalization condition with P = π, K and the chiral scale m 03 . The pion and kaon DAs have been determined at the scale 1 GeV in a recent global analysis [31] based on LO PQCD factorization formulas, which is at the same level of accuracy as this work. The results are quoted as where the Gegenbauer polynomials are defined as and the Gegenbauer moments in Eq. (16) are summarized as follows, a π 2 = 0.64 ± 0.08, a π 4 = −0.41 ± 0.10, a π 2P = 1.08 ± 0.15, a π 2T = −0.48 ± 0.33, a K 1 = 0.33 ± 0.08, a K 2 = 0.28 ± 0.10, a K 4 = −0.40 ± 0.07, a K 2P = 0.24, a K 2T = 0.35.
Note that the twist-3 DAs φ P K and φ T K , which were not obtained in Ref. [31], come from sum-rule calculations [49]. As stated before, we focus on the P -wave components in Eq. (1) proportional to P l=1 (2ζ − 1) = 2ζ − 1. The corresponding light-cone matrix element for a longitudinal meson pair is decomposed, up to the twist 3, into [30] where the two-meson DAs for the ππ, KK and Kπ pairs are parametrized as The Gegenbauer moments a 0,s,t 2ρ , a 0 1K * ,2K * ,4K * , and a 0 2φ will be determined in a global analysis in the next section. Since the current data are not yet precise enough for fixing the Gegenbauer moments in the twist-3 DAs φ s,t Kπ and φ s,t KK , they have been set to the asymptotic forms.
The elastic rescattering effects in a final-state meson pair can be absorbed into the time-like form factors F ,⊥ (ω 2 ), namely, the leading Gegenbauer moment B 00 (ω 2 ) in a two-meson DA according to the Watson theorem [50]. The resonant contribution from a ρ meson with a broad width is usually parameterized as the GS model [36] based on the Breit-Wigner (BW) function [51] in experimental investigations of three-body hadronic B meson decays, which interprets observed structures beyond the ρ(770) resonance in terms of heavier isovector vector mesons. Taking the ρ-ω interference and excited state contributions into account, we have the form factor [10,29,30] where m ρ,ω,j (Γ ρ,ω,j ), j = ρ ′ (1450), ρ ′′ (1700) and ρ ′′′ (2254), are the masses (decay widths) of the series of resonances, and c ω,j are the weights associated with the corresponding resonances. The function GS ρ (s, m ρ , Γ ρ ) is given by with the factors where the functions k(ω 2 ), h(ω 2 ) and β π (ω 2 ) are expressed as The function BW ω (ω 2 , m ω , Γ ω ) for the ω resonance takes the standard BW form [51]. We apply the RBW line shape for contributions from the intermediate resonances K * and φ of narrow widths to the form factors [7,8,11,35], with the mass-dependent widths where the masses m K * ,φ and the widths Γ K * ,φ of the K * and φ resonances, respectively, take the values in [37]. The magnitude of the spatial momentum of the meson P 1 , with the Källén function λ(a, b, c) = a 2 + b 2 + c 2 − 2(ab + ac + bc), is measured in the rest frame of the resonance, and | p 0 | is its value at the resonance mass. The orbital angular momentum L R in the two-meson system is set to L R = 1 for a P -wave state. Due to the limited knowledge on the form factors F ⊥ (ω 2 ), we assume the ratio [30], i = ρ, K * and φ, with f T i (f i ) being the tensor (vector) decay constants of the intermediate resonances.

A. Global Fit
We specify the parameters adopted in the numerical analysis below, including the masses (in units of GeV) [37] m B = 5.280, m Bs = 5.367, m b = 4.8, m K ± = 0.494, m K 0 = 0.498, m π ± = 0.140, m π 0 = 0.135, The Wolfenstein parameters in the Cabibbo-Kobayashi-Maskawa (CKM) matrix take the values in Ref. [52]: A = 0.836 ± 0.015, λ = 0.22453 ± 0.00044,ρ = 0.122 +0.018 −0.017 andη = 0.355 +0.012 −0.011 . We stress that ω B (s) in the B (s) meson DA, as an overall parameter, cannot be determined in our global analysis, but must be treated as an input. This is why we take its value extracted from the B (s) meson transition form factors in lattice QCD and light-cone sum rules, which has been also verified by the global study of two-body charmless hadronic B meson decays in [31]. The value of ω B (s) , together with the corresponding pion and kaon DAs fixed in [31], are then input into the present work on the three-body B meson decays for consistency. If the shape parameter ω B (s) is changed, the pion and kaon DAs need to be refitted accordingly, before they can be employed to constrain the two-meson DAs. Fortunately, the variation of ω B (s) causes less than 30% uncertainties for most of the branching ratios and negligible effects on the direct CP asymmetries as seen later, and is thus not expected to make a significant impact on the determination of the two-meson DAs. Hence, we focus only on the uncertainties of the Gegenbauer moments in the two-meson DAs propagated from experimental data here.
Equation (20) suggests that the total amplitudes A for the B (s) → P (ππ, πK, KK) decays with P = π, K can be expanded in terms of the Gegenbauer moments of the two-meson DAs. As a result, we decompose the squared amplitudes into the linear combinations of the Gegenbauer moments a 0,s,t 2ρ , a 0 1K * ,2K * ,4K * and a 0 2φ , and their products. We then compute the coefficients M , which involve only the Gegenbauer polynomials, to establish the database for the global fit.
Similar to the proposal in [31], we determine the Gegenbauer moments of the two-meson DAs by fitting the formulas in Eq. (30) with the Gegenbauer-moment-independent database to the measured branching ratios B and direct CP asymmetries A CP of the B (s) → P (ρ →)ππ, B (s) → P (K * →)Kπ and B (s) → P (φ →)KK decays. We adopt the standard nonlinear least-χ 2 (lsq) method [53], in which the χ 2 function is defined for n pieces of experimental data v i ± δv i with the errors δv i and the corresponding theoretical values v th i as The theoretical values v th i are the functions of the fitted Gegenbauer moments in the two-meson DAs. The lsq fit attempts to find the smallest χ 2 by adjusting the fitted parameters that bring the theoretical values closest to the data. The data with errors are treated as of the Gaussian type automatically in the fit packages, and the errors of the fitted parameters and of the theoretical values v th i come from experimental uncertainties by error propagation. To minimize statistical uncertainties, we should include maximal amount of data in the fit. On the other hand, those measurements with significance lower than 3σ do not impose stringent constraints, and need not be taken into account in principle. The data of those modes, which are affected by subleading contributions manifestly based on the previous PQCD studies [32,33], are also excluded, even though they may have higher precision. The B 0 → π 0 (ρ 0 →)ππ decay, dominated by the color-suppressed tree amplitude that is expected to receive substantial higher-order corrections [54], is a typical example.

B. Results
The Gegenbauer moments a 0 2ρ , a s 2ρ and a t 2ρ for the twist-2 and twist-3 ππ DAs in Table I are obtained from the fit to eight pieces of B → P (ρ →)ππ data marked by " †" in Tables II and III with χ 2 /d.o.f. = 2.6, whose errors mainly arise from experimental uncertainties. We point out that the measured B + → π + (ρ 0 →)ππ branching ratio, imposing a strong constraint on the Gegenbauer moment a s 2ρ , is considered in our fit, but the corresponding B + → π + ρ 0 data were excluded in the global analysis of two-body hadronic B meson decays [31]. It is seen that our Gegenbauer moments differ from the corresponding ones of the ρ(770) meson DAs derived in QCD sum rules [55] as mentioned before: the ππ DAs contain the ρ-ω mixing effect and the contributions from higher ρ resonances with finite widths via Eq. (21), so they need not to be the same as the ρ(770) meson DAs. Our Gegenbauer moments also differ from a 0 2ρ = 0.25, a s 2ρ = 0.75 and a t 2ρ = −0.60 chosen in Ref. [30] for two reasons at least. First, only the B → K(ρ →)ππ data were employed to constrain the ππ DAs in [30], while the additional B → π(ρ →)ππ data are included in our global analysis. Second, some B → K(ρ →)ππ data have been updated in this work.
A single Gegenbauer moment a 0 2φ is introduced into the KK twist-2 DA, and the twist-3 ones have been set to their asymptotic forms, since only two pieces of data from the B → K(φ →)KK decays in Table IV meet the required precision. The value of a 0 2φ , determined with χ 2 /d.o.f. = 0.35, is distinct from, but still consistent with that of the φ meson DA in QCD sum rules [55] within theoretical errors. Note that our a 0 2φ deviates from the value −0.50 ± 0.10 adopted in Ref. [4], where B s meson decays into charmonia plus a kaon pair were investigated. The deviation is understandable, because the choice of a 0 2φ depends on models for the uncertain charmonium DAs, as relevant data were accommodated.
The Kπ DAs are determined in a fit to six pieces of B (s) → P (K * →)Kπ data in Tables V and VI. We first work on Scenario I, in which the two Gegenbauer moments a 0 1K * and a 0 2K * of the twist-2 two-meson DA are fitted with χ 2 /d.o.f. = 1.5, and observe that a 0 2K * is slightly larger than unity as shown in Table I. A larger moment is not favored in view of the convergence of the Gegenbauer expansion. Therefore, one more Gegenbauer moment a 0 4K * is added in Scenario II, and a fit with The resultant a 0 2K * decreases a bit but with amplified uncertainty , and a 0 4K * is smaller than the unity. The Kπ branching ratios cannot give an effective constraint due to their larger experimental errors, such that the uncertainties of the Gegenbauer moments increase dramatically in Scenario II. For a similar reason, the obtained Gegenbauer moments differ from those of the K * meson DA in QCD sum rules [55], and from a 0 1K * = 0.2 and a 0 2K * = 0.5 chosen in the PQCD study on the B (s) → ψKπ decays [56]. We state that the fits based on the currently available data cannot discriminate the two scenarios effectively. As experimental precision increases, we will be able to impose more stringent constraints on those two-meson DAs.
The experimental data for comparison are quoted from Ref. [37]. Those data marked by † are included in the fit. The theoretical errors are attributed to the variations of the shape parameter ωB (s) in the B (s) meson DA and the decay constant fB (s) , of the Gegenbauer moments in the two-pion DAs, and of the hard scale t and the QCD scale ΛQCD.

Modes
Results With the fitted Gegenbauer moments in Table I, we calculate the CP averaged branching rations B and the direct CP asymmetries A CP in the LO PQCD formalism, and present the results in the central columns of Tables II-VI. The first theoretical uncertainty originates from the shape parameter ω B = 0.40 GeV or ω Bs = 0.48 GeV with 10% variation, and the decay constant f B (s) . The second one is from the Gegenbauer moments in the two-meson DAs. The last one is caused by the variations of the hard scale t from 0.75t to 1.25t, which characterizes the effect of next-to-leading-order QCD corrections, and of the QCD scale Λ QCD = 0.25 ± 0.05 GeV. The errors attributed to the CKM matrix elements are tiny and can be ignored safely. Note that the data for the B 0 → π + (ρ − →)ππ and B 0 → π − (ρ + →)ππ branching ratios in Table III represent the sum over these two modes. It is also the case for the measured B 0 s → K + (K * − →)Kπ and B + → π 0 (K * + →)Kπ branching ratios, and for the measured B 0 → π 0 (K * 0 →)Kπ and B 0 s →K 0 (K * 0 →)Kπ branching ratios in Table V.
One can also assess the uncertainties from the Gegenbauer moments a π 2,4 , a π 2P (T ) and a K (1,2,4) in the pion and kaon DAs. Taking the B + → π + (ρ 0 →)ππ and B 0 s → K − (K * + →)Kπ branching rations in Scenario II as examples, we obtain, given the errors in Eq. (18), It is seen that the former (latter) is more sensitive to the variation of the moment a π 4 (a K 2 ) in the twist-2 pion (kaon) DA. We remark that the total errors, derived by adding the individual ones from the moments in the pion and kaon DAs in quadrature and associated with the labels a π and a K below, respectively, are minor (less than 5%) compared with other uncertainties listed in Tables III and V: Therefore, the variation of the Gegenbauer moments in the pion and kaon DAs has little impact on the determination of the two-meson DAs. It is found that most of the considered data in Tables II and III are well reproduced, in particular those with higher precision. Larger deviation from the data is observed in the B + → π + (ρ 0 →)ππ and B + → π 0 (ρ + →)ππ branching ratios. It is ascribed to the involved color-suppressed tree contributions, which receive sizable next-to-leading-order corrections. The observables removed from the fit are also predicted in the LO PQCD formalism, and compared with the data in Tables II and III. Our prediction for the B 0 → π 0 (ρ 0 →)ππ branching ratio, which suffers significant subleading corrections as stated before, is still below the data, similar to that derived in the framework for two-body decays. Most of the A CP data for the B (s) → P (ρ →)ππ decays with P = π, K are not yet precise enough. We mention that A CP in the B + → π + ρ 0 mode has been predicted to be large and negative in most QCD approaches [10,57], including the present analysis on three-body decays as shown in Table III. However, its data are as small as 0.009 ± 0.019 [37]. Both the theoretical and experimental errors need to be reduced greatly in order to tell whether the discrepancy really stands as a puzzle.
Both the B → K(φ →)KK data considered in the fit are well reproduced with a single Gegenbauer moment a 0 2φ as indicated in Table IV. Our predictions for the branching ratios and direct CP asymmetries excluded in the fit, mainly associated with B s meson decays, can be confronted by more precise data in the future. All the available A CP data for the B → P (φ →)KK decays with P = π, K have large errors. The central values of the prediction and the data for the B + → π + (φ →)KK branching ratio are different, but still agree with each other within uncertainties.  Table II but for the B (s) → P (φ →)KK decays with P = π, K.

Modes
Results Overall speaking, Scenario II reproduces the considered B (s) → P (K * →)Kπ data with P = π, K slightly better than Scenario I does as seen in Tables V and VI. The B s → P (K * →)Kπ branching ratios differ between the two scenarios more than the B → P (K * →)Kπ branching ratios do. This feature is understandable, because the former involve the B s → (K * →)Kπ transition form factors, which are more sensitive to the variation of the Gegenbauer moments in the Kπ DA.
Hence, more precise B s → P (K * →)Kπ data are crucial for fixing the Kπ DAs. The direct CP asymmetries A CP in some B (s) → P (K * →)Kπ modes depend on the chosen scenarios strongly, implying that more accurate Kπ DAs are necessary for predicting these observables unambiguously. The central value of the predicted B 0 s → π + (K * − →)Kπ branching ratio in Scenario II, which is already much lower than in Scenario I, remains above the data. It deserves more thorough theoretical and experimental investigations. Similarly, most of the A CP data for the B (s) → P (K * →)Kπ decays have substantial uncertainties so far, so it is not yet possible to make a meaningful comparison with our results.   16) and (20), respectively, are expanded up to the fourth-order Gegenbauer polynomial without the third-order term, which exists in general. We find that the SU (3) breaking effects in the considered decays can be well accounted for by the first-order term alone. That is, even the third-order term is included into the fit, its value turns out to be small, and does not modify the fit outcomes much. Taking the eight B (s) → P (K * →)Kπ (P = π, K) decays as examples, we perform the global fit with a 0 3K * being included, and obtain the Gegenbauer moments of the Kπ twist-2 DA and the branching ratios in Table VII. It is seen that the central value of a 0 1K * increases only a bit with an enlarged uncertainty and a 0 2K * stays the same, compared with those from Scenario I in Table I, and the central value of a 0 3K * is tiny. The corresponding branching ratios in Table VII also change very little, compared with those in Tables V and VI. The above observations support that the SU (3) breaking effects in the considered modes can be explained by the a 0 1K * term alone under the current data precision. Hence, the neglect of the a 0 3K * term in this work is justified. Besides, it is not practical to include many parameters in the fit because of the limited amount of experimental data at present. For a similar reason, the asymptotic forms of the Kπ twist-3 DAs φ s,t Kπ are adopted in our analysis. The same argument applies to the expansion of the kaon DAs in Eq. (16), where the higher moments responsible for SU (3) symmetry breaking effects are also absent. We will explore the impact of these neglected Gegenbauer polynomials systematically in the future, when more experimental data with improved precision are available.  [37]. For simplicity, only the theoretical errors from the Gegenbauer moments are presented.

Modes
Results It is noticed that the parametrization of the parton momenta in Eqs. (2) and (3) introduces the dependence on the light meson mass m 3 into the hard kernels and the Sudakov exponents, as explicitly shown in the Appendix. Since both these factors are perturbative pieces in a PQCD factorization formula, they should be insensitive to a light scale. Therefore, we test the sensitivity of our numerical results to this light scale by setting it to zero in the hard kernels and the Sudakov exponents. The corresponding branching ratios and direct CP asymmetries for two typical modes, B + → K + (ρ 0 →)ππ and B 0 → π − (K * + →)Kπ in Scenario II, are presented in Table VIII. The neglect of the kaon mass for the former mode causes about 10% variation in the branching ratio and the direct CP asymmetry. The quantities associated with the latter mode are relatively stable with respect to the neglect of the pion mass as expected. The insensitivity to the light scale confirms that our parametrization for kinematic variables in three-body B meson decays is reasonable. VIII: CP averaged branching ratios and direct CP asymmetries of the B + → K + (ρ 0 →)ππ decay and the B 0 → π − (K * + →)Kπ decay in Scenario II with and without the light meson masses in the hard kernels and the Sudakov exponents. The experimental data are quoted from [37]. The sources of the theoretical errors are the same as in Table II.

Modes
Results (with light mass) Results (without light mass) Data We make a more aggressive attempt in this subsection to determine the dependence of the Gegenbauer moments in the twomeson DAs on the meson pair invariant mass. As stated in the Introduction, it is unlikely to extract the exact dependence from current data, so we simply expand the Gegenbauer moments up to the first power in ω 2 , and examine whether the additional linear terms can be constrained effectively in the global fit. Consider the parametrizations of the di-pion DAs, with the free parameters a 0,s,t 2ρ and c 0,s,t ρ . The above parametrization follows the power series for the ω 2 -dependent Gegenbauer moments derived in Ref. [22].
The global fit to the same set of B (s) → P (ρ →)ππ data with P = π, K leads to the outcomes in Table IX with a smaller χ 2 /d.o.f. = 0.51, which are not difficult to understand: varying ω 2 around the ρ resonance in its width window, we find that the values of a 0,s,t 2ρ (1 + c 0,s,t ρ ω 2 ) are in fact consistent with the corresponding ones in Table I. The consistency is particularly obvious for a t 2ρ (1 + c t ρ ω 2 ) with the tiny coefficient c t ρ . It is observed from Table IX that the parameters for the twist-3 DA φ s ππ , which gives sizable contributions to branching ratios, can be constrained effectively by the current data. It suggests that the determination of the ω 2 -dependent Gegenbauer moments is promising, when more precise data are available in the future. Because our purpose is to demonstrate the potential to extract the ω 2 dependence of the Gegenbauer moments, we will not work on the Kπ and KK DAs. The effect of including the ω 2 dependence of the Gegenbauer moments is similar to that of introducing more parameters. That is, the fit quality is improved with a lower χ 2 /d.o.f. at the cost of larger uncertainties for fit results as shown in Table X. For example, the reproduced branching ratios for the B + → K + (ρ 0 →)ππ and B + → π + (ρ 0 →)ππ decays get closer to the data, which have relatively higher precision. However, the uncertainty caused by the variation of the di-pion DAs is amplified compared to the second source of errors in Table II.

IV. CONCLUSION
In this work we have performed a global fit of the Gegenbauer moments in two-meson DAs to measured branching ratios and direct CP asymmetries in the three-body hadronic B meson decays B → V P 3 → P 1 P 2 P 3 with V = ρ, φ, K * and P 3 = π, K in the LO PQCD approach. Two-meson DAs, collecting both nonresonant and multi-resonance contributions, serve as crucial nonperturbative ingredients of factorization theorems for the above decays. The Gegenbauer moments of the pion and kaon DAs determined in the LO global analysis of two-body hadronic B meson decays have been input for theoretical consistency. To facilitate the numerical study, we have constructed a Gegenbauer-moment-independent database, via which a decay amplitude is decomposed into a linear combination of the relevant Gegenbauer moments in the two-meson DAs. It was noticed that the fitted Gegenbauer moments differ from those associated with an intermediate resonance which decays into the meson pair, and from those adopted in previous PQCD calculations. This observation indicates that the Gegenbauer moments of a two-meson DA cannot be inferred from sum-rule results for an intermediate resonance, and their global determination is essential.
We have examined two scenarios for the determination of the Kπ DAs in order to check the convergence of the Gegenbauer expansion, and the sensitivity of the fitted observables to our setup. It was found that the Gegenbauer expansion is improved by increasing the number of Gegenbauer moments at the cost of large uncertainties for fit outcomes, and that the branching ratios of B s meson decays and direct CP asymmetries in some modes are more sensitive to the chosen scenarios. Hence, more accurate Kπ DAs are necessary for predicting these quantities unambiguously. We state that our fits have not been able to discriminate the two scenarios effectively. We have also explored the potential to fix the dependence of the Gegenbauer moments on the meson-pair invariant mass, and confirmed that at least the parameter for the twist-3 DA φ s ππ can be constrained to some extent by the current data. Therefore, the determination of the dependence on the meson-pair invariant mass is promising, when data become more precise.
We mention that the three-body charmless hadronic B meson decays included in this work have been studied in Refs. [6,7,10,11] in a scattered manner. The improvements compared to the earlier studies contain (1) the partonic kinematic variables have been refined to take into account finite masses of final-state mesons, such that the SU (3) symmetry breaking effects in the decays can be evaluated more precisely; and (2) the Gegenbauer moments in the two-meson DAs have been determined in a global analysis for the first time, which are valuable for future applications of the PQCD framework to multi-body B meson decays; (3) the dependence of the Gegenbauer moments on the meson-pair invariant mass has been probed for the first time. Because of (1), the numerous hard kernels involved in the various modes need to be modified, which have been presented, together with the factorization formulas for the decay amplitudes, in the Appendix. The refined partonic kinematics is general enough for its extension to multi-body B meson decays into arbitrary massive final states. For (2), we remind that different Gegenbauer moments for the Kπ DAs were taken in the previous scattered studies, such as Refs. [7] and [11], so our work facilitates a consistent understanding of multi-body B meson decays. We have shown that the preferred central value of, for instance, the Gegenbauer moment a 0 1K * is 0.31, instead of 0.2 in [7] or 0.05 in [11] (but note the large theoretical uncertainties).
It has been demonstrated that most of the data considered in the fit are well reproduced, namely, the fit quality is satisfactory. It implies that the two-meson DAs presented in this paper are ready for applications to other multi-body hadronic B meson decays involving the same meson pairs. With the obtained Gegenbauer moments, we have made predictions for those observables, whose data were excluded in the fit because of their substantial experimental errors or significant subleading contributions to the corresponding factorization formulas. Except the B 0 s → π + (K * − →)Kπ branching ratio, our predictions agree with the data within uncertainties in the former case. Since our results were still derived in the LO PQCD approach, the data in the latter case remain unexplained, and deserve more through investigations. As pointed out before, the precision of the extracted two-meson DAs can be improved systematically, when higher-order and/or higher-power corrections to three-body hadronic B meson decays are taken into account in our formalism. At the same time, more precise measurements are urged, especially those of CP asymmetries. These efforts will strengthen the constraint on the Gegenbauer moments and sharpen the confrontation between theoretical predictions and experimental data.