Extraction of the CKM phase $\gamma$ using charmless 3-body decays of $B$ mesons

The weak phase $\gamma$ is extracted from three-body charmless decays of $B$ mesons following a method proposed by Bhattacharya, Imbeault \&London. The result is obtained by combining the BABAR amplitude analyses for the processes $B^0 \rightarrow K^+ \pi^0 \pi^-$, $B^0 \rightarrow K_S^0 \pi^+ \pi^-$, $B^0 \rightarrow K_S^0 K_S^0 K_S^0$, $B^0 \rightarrow K^+ K_S^0 K^-$ and $B^+ \rightarrow K^+ \pi^+ \pi^-$, under the assumption of $\alpha_{\rm{SU(3)}}$ flavour symmetry. Six possible solutions are found. One solution is compatible with the Standard Model while the others are not. It is also found that, when averaged over the entire Dalitz plane, the effect of $\alpha_{\rm{SU(3)}}$ breaking on the analysis is only at the percent level.


Introduction 2 Method of extraction of γ
We begin with a review of the method for extracting γ from B → Kππ and B → KKK decays. There are several ingredients. First, the amplitudes for three-body B → P P P decays can be written in terms of diagrams [6,10]. These diagrams are similar to those of two-body decays [11,12], except that here it is necessary to "pop" a quark pair from the vacuum. Also, in contrast to two-body diagrams, the three-body diagrams are momentum dependent.
Second, one can fix the symmetry of the final state in B → P 1 P 2 P 3 by using its Dalitz plot [6]. We define the three Mandelstam variables s i ≡ (p j + p k ) 2 , where p i is the momentum of P i , and ijk = 123 231 or 312. (These obey s 1 + s 2 + s 3 = m 2 B + m 2 1 + m 2 2 + m 2 3 .) Experimentally, one can reconstruct the decay amplitude A (s 1 , s 2 ), which varies as a function of position in the Dalitz plot. The amplitude that is fully symmetric under permutations of the final-state particles is then given by The symmetrised amplitude A fs has a sixfold symmetry in the Dalitz plane. In effect, the plane can be divided into six regions; the structure and information in each region is identical to the others. It is therefore sufficient to consider points in one sixth of the symmetrised Dalitz plane. Third, in Ref. [7] it was shown that, as is the case in two-body decays [13,14,15], under flavour SU (3) there are relations between the electroweak penguin (EWP) and tree diagrams for b → s transitions. For the fully-symmetric final state, these take the form where the c i are Wilson coefficients and λ (s) p ≡ V * pb V ps (the V ij are elements of the CKM matrix). The method uses B → Kππ and B → KKK decays. In the B → Kππ diagrams, the quark pair popped from the vacuum is uū or dd (under isospin, these diagrams are equal). On the other hand, the B → KKK diagrams have a popped ss pair. Now, under flavour SU(3) symmetry, which is required for the EWP-tree relations eq. (2), diagrams with a popped ss quark pair are equal to those with a popped uū or dd. In other words, under SU(3) the diagrams in B → KKK decays are the same as those in B → Kππ decays.
Of course, flavour SU(3) symmetry is not exact, so one must keep track of SU(3) breaking. Technically, there will be an SU(3)-breaking factor for each diagram. However, if all such quantities are included, there will be too many unknown parameters to perform a fit. For this reason, we make the assumption that the size of SU(3) breaking is the same for all diagrams, so there is a single SU(3)-breaking parameter α SU(3) relating B → Kππ and B → KKK decays (α SU(3) = 1 corresponds to the flavour-SU(3) limit). The idea behind this assumption is as follows. As noted above, the diagrams are momentum dependent. This means that the size of SU(3) breaking associated with a particular diagram, α D , varies from point to point on the Dalitz plot. Specifically, α D will be > 1 at some points and < 1 at others. When one averages over all points, α D − 1 will be small, and this will be true for all diagrams. For this reason, we make the assumption that the size of SU(3) breaking is the same for all diagrams, and we expect it to be small. As we will see, in our fits α SU(3) − 1 is found to be at the percent level, which supports our assumption.
Three B → Kππ and two B → KKK decays are used in this analysis. They are As shown in Ref. [9], when the EWP-tree relations of eq. (2) are used, the fully-symmetric amplitudes for the five modes can be expressed as linear combinations of five effective diagrams: Here the complex parameters A, B, C, D andP uc are linear combinations of momentum-dependent diagrams (and will, in general, vary across the phase space), γ is the CKM angle of interest, and κ is a constant defined in eq. (2). As noted above, the real quantity α SU(3) parametrizes the breaking of flavour SU(3) symmetry and can also vary across the Dalitz plot. In the absence of any SU(3)-breaking effects, α SU(3) = 1, so that the amplitudes of B + → K + π + π − and B 0 → K + K 0 K − modes are equal. In this limit, the fifth decay mode provides no additional information and can be dropped from the analysis. For each decay mode, a set of linearly-independent observables can be formed: whereĀ fs (s 1 , s 2 ) denotes the fully-symmetric amplitude of the conjugate process. The observables X, Y , and Z are related to the effective CP -averaged branching fraction, the direct CP asymmetry, and the indirect CP asymmetry. For a given decay, their values depend on the position in the Dalitz plane. The observable Z has no physical meaning for flavour-specific final states such as K + π 0 π − and K + π + π − . In this study, we take as experimental inputs the amplitude models obtained by BaBar in Refs. [16,17,18,19,20]. The BaBar analysis of B 0 → K 0 S K 0 S K 0 S [19] was time-integrated and CP -averaged; since no distinction was made between B 0 and B 0 , only the observable X is accessible for this mode. As noted in Ref. [9], this implies a simplification in the expression for its amplitude compared with eq. (3). To be specific, the requirement that Y = Z = 0 implies thatP uc = 0, so that Since for each mode the observables X, Y, Z depend upon the fully-symmetric amplitude A fs [Eq. (4)], and A fs is related to the theory parameters by Eqs. (3) and (5), the observables may be written as functions of those theoretical parameters. Expressing them in terms of magnitudes and strong phases (U = ue iφu for U = A, B, C, D), and setting φ a = 0 without loss of generality, the following relations are obtained: If γ is extracted at a single point (s 1 , s 2 ) on the Dalitz plane, there are nine real, unknown parameters: four magnitudes (a, b, c, d), three strong phases (φ b , φ c , φ d ), γ, and α SU (3) . From the experimental input, there are eleven observables: three (X, Y, Z) for each of the modes K 0 S K + K − and K 0 S π + π − , two (X, Y ) for each of the modes K + π + π − and K + π + π 0 , and one (X) for K 0 is fixed to unity, there are instead eight unknown parameters and nine observables. In both cases, there are more observables than theory parameters, and γ may be extracted with a fit.
One can instead determine γ using information from several points on the Dalitz plane simultaneously. This increases the number of observables, but it also increases the number of unknowns, since all those parameters that describe the strong dynamics of the decay can vary across the Dalitz plane (whereas γ of course does not). Thus, if N points on the Dalitz plane are used, there are 11N observables and 8N + 1 unknown parameters when α SU(3) is allowed to vary, or 9N observables and 7N + 1 unknowns if α SU(3) is fixed to unity. For any N ≥ 1, the number of observables exceeds the number of unknowns in both cases, so that γ may again be extracted with a fit.

Practical implementation of the extraction method
The five BaBar analyses use the isobar formalism to parametrise the variation of the amplitude across the Dalitz plane. In this approach, the total amplitude at a point (s 1 , s 2 ) in the Dalitz plane is given by the coherent sum of the amplitudes of n individual decay channels: where the isobar coefficients, c j , are complex numbers containing all the weak phase dependence, and the lineshapes, F j , are wave functions (such as Breit-Wigner functions) that describe the dynamics of the decay amplitudes. The isobar coefficients and the lineshapes given in BaBar's papers are used to compute the amplitudes of the different decay modes as a function of position in the Dalitz plane. This calculation is implemented with the LAURA ++ software package [21]. The uncertainties on the isobar coefficients quoted by BaBar, and the associated correlation matrices, are used to compute the experimental uncertainties on the extracted values of the angle γ. For the decay mode B 0 → K 0 S K 0 S K 0 S , no correlation matrix was quoted, and the correlations are therefore neglected.
After encoding the isobar model for a decay mode in LAURA ++ , the amplitudes for the process and its CP conjugate may be calculated at any point, or set of points, on the Dalitz plot. From the amplitudes, those observables (X, Y , Z) that are well-defined may be computed according to eq. (4). The uncertainties on the isobar coefficients are propagated to obtain the uncertainties on the observables along with their correlations.
It is also possible to compute the observables from the theory parameters eq. (6); for a set of N points, those theory parameters consist of γ plus N instances of the amplitude and SU(3)-breaking parameters, The χ 2 may then be computed for compatibility between the observables expected given these theory parameters, and the observables obtained from experimental inputs. To constrain γ, a χ 2 scan is performed as follows: γ is fixed to a certain value, then a χ 2 minimisation is performed with the other 8N parameters free to vary. In principle, for a global minimisation the final values of the fit parameters and the χ 2 should not depend on the initial values of the parameters. In practice, for multidimensional fits some dependency is observed due to the presence of secondary local minima. In order to obtain a more robust estimate of the global minimum the minimisation process is repeated 500 times, with values of the initial parameters varied randomly in a large physical range each time. Then, for each fixed value of γ, the smallest χ 2 is retained. The value of γ is then increased by one step and the minimisation repeated. Performing this many times, a scan of the χ 2 as a function of γ is obtained. The minima of this scan are the preferred values for γ. The procedure for finding the minima in a scan is detailed in Appendix A. The asymmetric statistical uncertainty on each solution is then estimated as the change in γ required to produce a change of one unit in χ 2 from the minimum.
Flavour SU(3) symmetry can be broken locally when considering single points on the Dalitz plane. However, as shown in Sec. 6 below, SU(3)-breaking effects are small when averaging over a large number of points. For this reason, as well as to minimise the statistical uncertainties and to use the maximum amount of information possible, it is desirable to extract γ using the largest possible number of points. In principle, γ could be obtained from an arbitrarily large set of points. However, the observables can be highly correlated between points, especially if the same resonance (or non-resonant component) is the dominant contributor to the points in one of the decay modes. High correlations have an impact on the covariance matrix, which becomes approximately singular and not invertible. This imposes practical limitations on the choice of points: the total number of points that can be used in a fit is finite and small. For these modes and isobar models, it is found that no more than three points can be used, and not all three-point combinations are possible. In order to avoid dependence on the choice of points (i.e., to avoid experimenter's bias), the scan procedure is carried out repeatedly with random combinations of three points, applying a filter to reject sets of points for which the correlation matrix contains entries above 70%. In total, 501 random sets of three points that pass this filter are used.
For each of the solutions (scan minima) for γ, the final result is taken to be the average over the central values for that solution, and the uncertainty on that result is the average of the uncertainties from individual scans. Note that, fluctuations aside, the average uncertainty does not decrease as more scans are added.
As discussed in Sec. 2, the extraction of γ can be performed with four modes (fixing α SU(3) = 1 and so neglecting flavour SU(3) symmetry-breaking effects), or with five modes (allowing α SU(3) to vary and parametrising those effects). For reasons of fit stability and convergence, the method with four modes is chosen to be the baseline for the results and will be presented in the next section. The method using five modes is then used to assess the systematic uncertainties associated with SU(3) breaking; the corresponding uncertainties are given in Sec. 5, with the results for the five-mode method described in more detail in Appendix B.
4 Results with four modes, α SU (3) = 1 The method described in the previous section is applied, and χ 2 scans for γ are obtained. Six distinct minima are found; averaging over the 501 sets of three points in the Dalitz plot, their central values µ and asymmetric experimental uncertainties (σ L , σ R ) are given in Table I. The experimental uncertainties are below 11 • in each case. The third of these minima, at 68.9 • , is compatible with the current world-average value of γ [2].  Since each combination of three points carries different information, the form of the χ 2 scans vary from one combination of points to the next. The central values fluctuate, and, in some instances, not all of the six minima are present. The distribution of the minima across the scans is shown in Figure 1, and the rates at which the minima are found are given in Table II; each of them is found in more than 90% of scans.

Systematic uncertainties
The experimental statistical and systematic uncertainties on the amplitude models used as inputs are already included in the results given in Table I. Two additional sources of systematic uncertainty, discussed below, are considered in this study. The first relates to the combination of the minima obtained with different sets of three points in the Dalitz plot. The second relates to flavour SU(3) breaking. The results are summarised in table III. The form of the χ 2 scan varies according to the points chosen, and in some instances a minimum is found successfully but is not well separated from another nearby minimum, such that if the two minima are at µ 1 and µ 2 with χ 2 values χ 2 1 and χ 2 2 , no value of γ in the range µ 1 < γ < µ 2 (or µ 2 < γ < µ 1 ) has a χ 2 value greater than or equal to χ 2 1 + 1. This means that the algorithm set out in Sec. 3 cannot determine the experimental uncertainty on µ 1 . These minima are referred to as poorly resolved and are not included in the average from which the overall results are obtained (Table I). Discarding these minima could have a systematic effect on the average (e.g., if µ 1 < µ 2 are two nearby minima then upward fluctuations in µ 1 are more likely to be too close to µ 2 to resolve than downward fluctuations in µ 1 , potentially causing a negative bias in µ 1 ). To assess this effect, the analysis is repeated including all minima from all scans in the average, even those that are not well resolved. The systematic uncertainty is then assessed as where µ is the central value obtained including only well-resolved minima in the average, and µ all is the central value obtained when including both well-resolved and not-well-resolved minima. The values obtained are given in Table III and are below 1.5 • for each minimum. The extraction performed with four modes does not take into account flavour SU(3) breaking. While it is not practical to allow for SU(3) breaking in a completely general way in this analysis, the scale of the effect can be assessed by allowing the SU(3)-breaking parameter α SU(3) to vary and seeing how much the values of γ change. To this end, the analysis is repeated using five modes instead of four, and with α SU(3) free to vary as an additional real parameter in the fit. As before, a χ 2 scan for γ is obtained with hundreds of random combinations of three points in the Dalitz plot, and for each scan the minima are found. (More details are given in Appendix B.) For each minimum, the central value of γ is averaged over the scans as before. These estimates using five modes (µ 5 modes ) may then be compared to the value for that minimum obtained with the baseline, four-mode procedure (µ) to assess how large an effect flavour SU(3) breaking has on the value of γ: The values obtained are given in Table III and are below 3 • for each minimum. More tests of the validity of the flavour SU(3) symmetry hypothesis are described in Sec. 6.

Studies of flavour SU(3) breaking
The assumption of flavour SU (3), and specifically that α SU(3) = 1, is tested in two further ways. The first involves comparing the amplitudes of two modes related by flavour SU(3) as a function of position in the Dalitz plane. The second consists of determining the value of α SU(3) over the Dalitz plane from fits to the amplitude models.
6.1 Comparison of the amplitudes of B 0 → K S K + K − and B + → K + π + π − From inspection of the last two lines of eq. 3, there is a linear relationship between the fully symmetric amplitudes for B 0 → K 0 S K + K − and B + → K + π + π − : The value of the parameter α SU(3) , a measure of the amount of the local flavour SU(3) breaking, can be inferred by comparing the values of the amplitudes of these modes at different points on the Dalitz plane [10]. We define the following ratio: R(s 13 , s 23 ) = A fs (B + → K + π + π − ; s 13 , s 23 ) + A fs (B − → K − π − π + ; s 13 , s 23 ) where A fs (X; s 13 , s 23 ) is the symmetrised amplitude for the decay mode X measured at point (s 13 , s 23 ). The ratio is an estimate of α SU(3) at that point. Figure 2 (a) shows the value of R(s 13 , s 23 ) as a function of position in the Dalitz plane. Significant deviations from unity are seen, especially near resonances. This is unsurprising, given that flavour SU (3) is broken by the mass difference between s and u/d quarks. A histogram of the values of R, sampled uniformly across the Dalitz plane, is shown in Fig. 2 (b). The distribution peaks near one, and the average value is 1.028, rather close to unity. This suggests qualitatively that, while SU(3) is strongly violated locally, it holds reasonably well when averaging across the phase space.

Fitted value of α SU(3) over the Dalitz Plane
Another approach is to determine α SU(3) from a fit. For this exercise, individual points in the Dalitz plane are considered (as opposed to sets of three points). A uniform grid of 386 points is used. For each point, a similar procedure is followed to that described in Sec. 3, with a χ 2 minimisation carried out with γ being fixed to a certain value γ i and the other physics parameters, including α SU(3) , being free to vary. As before, the fit is repeated 500 times (for each point) with the initial parameter values randomised, and the solution with the smallest χ 2 after the fit is retained. However, instead of scanning for γ across the full range, the exercise is only performed for six values γ i corresponding approximately to the six minima given in Table I. At each point and for each value of γ i tested, the value of α SU(3) for the best-fit solution is recorded.
Averaging over the uniform grid of points, the mean values of α SU(3) are given in Table IV for each γ i . Each value is close to unity, and negligible variation in the average α SU(3) is seen between the six minima. The variation of α SU (3) with position in the Dalitz plot is illustrated in Fig. 3, in which at each point in the symmetrised Dalitz plot the fitted values of α SU(3) from the six γ i are averaged. Similar structure is seen to that observed in Fig. 2, and the breaking of flavour SU(3) is clearly seen near resonances.

Conclusion
The method of extracting the weak phase γ from three-body charmless decays of the B meson developed by Bhattacharya, Imbeault and London [9] is applied to amplitude models of five charmless three-body decays of B mesons obtained by the BaBar collaboration [16,17,18,19,20]. Six solutions for γ are found:  (3) breaking as well as the impact of poorly resolved minima on the procedure. The statistical uncertainty is dominant, and is below 11 • for each of the six solutions. This is approximately a factor two larger than the uncertainty on the world-average value of γ, and allows the value obtained from these loop-level processes to be compared to the tree-dominated average. The presence of multiple solutions may reflect trigonometric ambiguities in the amplitudes. Further tests of the flavour SU(3) symmetry hypothesis were performed, studying the variation in the SU(3)-breaking parameter α SU(3) across the phase space. Strong local variation is seen, comparable to the ∼ 30% level typically considered, but the average value of α SU(3) is found to be close to 1 (corresponding to SU(3) symmetry) within a few percent.
The study presented in this paper is a complete proof of principle, including fully-propagated experimental uncertainties. It would benefit from additional and more precise experimental inputs; results from Belle II and LHCb would be welcome. It is worth noting that certain modes are well suited to the LHCb detector (e.g. B + → K + π + π − ), while others are better adapted to Belle II (e.g. B 0 → K 0 S K 0 S K 0 S ). Given this, one interesting possibility would be a simultaneous fit of the physics parameters to datasets of both experiments using a framework such as JFIT [22].
Further developments on the theoretical side would also be welcome, such as considering other symmetry states (fully antisymmetric or of mixed symmetry). This would add information, thereby reducing the statistical uncertainties, and might help to resolve the ambiguities and determine whether the value of γ found using loop-level processes is or is not equal to that obtained using tree-level decays.

Appendices A Algorithm for extracting the minima
The following algorithm is used to find the minima in a given scan: 1. Start at the first point.
2. Define the current window to be the range of γ spanned by the current point plus the next 19 consecutive points. Fit those 20 points with a 3 rd -order polynomial function.
3. Determine the minimum of the fitted polynomial (at x = γ min , y = χ 2 (γ min ) 4. Reject the minimum (γ min , χ 2 (γ min )) if any of the following is true: • The value of γ min is outside the window.
• The polynomial fit is of poor quality (its fit χ 2 is greater than 5).
5. Move along one point, then go back to step 2 (unless the points have been exhausted).
Usually, when a minimum is identified, it will be found by several consecutive polynomial fits (steps 2-4). Due to statistical fluctuations, the value of γ min will differ slightly between these; the average value is taken.

B Extraction of γ with five modes varying α SU(3) in the fit
The analysis (described in Sec. 3) was carried out using four decay modes and with α SU(3) fixed to unity as a baseline. To assess the systematic effect of SU(3) breaking, a similar procedure was used with five decay modes and with α SU(3) free to vary in the fit. The following changes were made to the procedure: the rejection criterion on the correlation between sets of points was relaxed from 70% to 80%, and the number of random set of three points was reduced from 501 to 401. The fit behaviour was found to be less stable, with convergence of the χ 2 minimisation in around 80% of cases (rather than 100% in the baseline). The frequency with which the minima were identified was also reduced (as shown in Table V). The reduced stability is taken to be due to the increased number of free parameters, and the consequent increase in the size of the covariance matrix.
The results of the procedure with five modes are shown in Table VI, giving the central values (µ), asymmetric experimental uncertainties (σ L , σ R ), and the recomputed systematic uncertainty due to poorly resolved minima (|µ − µ all |). The systematic uncertainty associated with SU(3) breaking is also given in the table; this is the same as before by construction. The distribution of the minima across the scans is shown in Figure 4. The results for the minima are compatible with the ones obtained with four modes.   (3) to vary in the fit. For each minimum, the central value for γ is given (µ), along with the asymmetric experimental uncertainty on the left-and right-hand sides (σ L , σ R ). The quantities |µ − µ all | and |µ 4 modes − µ 5 modes | are taken as estimates of the systematic uncertainties due to poorly resolved minima and flavour SU (3)