A novel unbinned model-independent method to measure the CKM angle $\gamma$ in $B^{\pm} \to DK^{\pm}$ decays with optimised precision

We present a novel unbinned method to combine $B^{\pm} \to DK^{\pm}$ and charm threshold data for the amplitude-model unbiased measurement of the CKM angle gamma in cases where the D meson decays to a three-body final state. The new unbinned approach avoids any kind of integration over the D Dalitz plot, to make optimal use the available information. We verify the method with simulated signal data where the D decays to $K_S \pi^+ \pi^-$. Using realistic sample sizes, we find that the new method reaches the statistical precision on gamma of an unbinned model-dependent fit, i.e. as good as possible and better than the widely used model-independent binned approach, without suffering from biases induced by a mis-modeled D decay amplitude.


Introduction
The precision measurement of the CP violating phase γ (ϕ 3 ) is one of the primary goals of flavour physics today. The measurement of γ in B ± → DK ± and related decays [1][2][3][4][5][6] has negligible theoretical uncertainty [7]. The current precision on γ is dominated by LHCb, which measures γ = 63.8 +3.5 −3.7 • [8,9]. The vast, clean datasets expected from the recently commissioned LHCb upgrade [10] and Belle II [11] will allow a sub 1 • precision on γ, LHCb upgrade II is expected to reduce this further to 0.35 • [12,13]. The parameter γ is therefore set to become the most precisely measured parameter in the Cabibbo-Kobayashi-Maskawa (CKM) description of CP violation in the quark sector [14,15], giving it a pivotal role in the search for new physics by over-constraining the Standard Model with precision measurements. To fully benefit from this potential, it is critical to control systematic uncertainties.
The measurement of γ in B − → DK − and its CP conjugate is possible because the neutral D meson in this decay is a superposition of D 0 and D 0 that depends on γ: D ∝ D 0 + r B e i(δ B −γ) D 0 , where r B ∼ 0.1 and δ B is a CP -conserving phase induced by the strong interaction. The neutral D meson is reconstructed in a final state accessible to both D 0 and D 0 . For multibody decays of the D, such as K 0 S π + π − , the CP -violating phase γ is obtained from analysing the amplitude structure of the D decay in both B − → DK − and B + → DK + transitions; this approach is known as the BPGGSZ method [3][4][5].
In order to measure γ in this way, the phase difference between the D 0 and D 0 decay amplitude across the Dalitz plot needs to be known. While this information can be obtained from amplitude models, these models have well-known shortcomings that make the phase information they provide unreliable, leading to significant, difficult-to-quantify systematic uncertainties. For this reason, model-independent methods are required to achieve the ultimate precision on γ. Similar considerations apply to measurements of D-mixing [16].
Model-independent methods currently in use rely on integrating over all or parts (bins) of the multibody phase space of the D decay [4,[17][18][19][20][21]. Two unbinned model-independent methods have been proposed recently. The method described in [22] is based on projecting the two-dimensional Dalitz plot down to one dimension, where the phase difference between the D 0 and D 0 amplitudes is parameterised as a Fourier series. The authors of [23] extract γ from set of cumulative functions defined across the Dalitz plot in a way that is inspired by the Kolmogorov-Smirnov test. This method is independent of the phase difference between the D 0 and D 0 amplitudes.
In all cases, one can expect some information loss due to the integration or projection process involved. Here we present a new unbinned model-independent method that optimally uses the full information across the two-dimensional Dalitz plot of a three-body decay of the neutral D meson. Using the example of D 0 → K 0 S π + π − , we will show that this method has, with current and plausible future data sample sizes, the potential to reach essentially the same statistical precision on γ as a model-dependent method, without suffering from the associated model uncertainty.
This paper is organised as follows: In Sec 2 we remind the reader of the formalism for the measurement of γ in B ± → DK ± decays, and use this opportunity to introduce our notation and phase convention. Section 3 describes the new quasi model-independent method introduced in this paper. In Sec 4, we evaluate the performance of the new method in simulation studies for the measurement of γ in B ± → DK ± , D → K 0 S π + π − decays. Finally, in Sec 5, we conclude.

Formalism
In this section we outline the formalism for the measurement of γ and the use of charm threshold data developed in [1][2][3][4][5][6], to remind the reader and fix the notation. For simplicity, we ignore the effects of D mixing [24][25][26][27]. In most cases, these effects are small enough to be negligible [26]. Where they are not, they can be taken into account as described in [27].

Notation and conventions
We use the notation A f D (p) for the decay amplitude of the D 0 meson to a final state f , at point p in the D → f Dalitz plot and A f D (p) for that of the CP -conjugate process.
Our phase convention for the CP operator is such that CP |D 0 ⟩ = |D 0 ⟩, which is the usual practice in the context of beauty decays to charm. An alternative convention where CP |D 0 ⟩ = −|D 0 ⟩ is widely used in the context of charm physics.
With our phase convention, and assuming the absence of CP violation in charm de- . For self-conjugate final states such as K 0 S π + π − , which we will use for the simulation studies presented below, p and p are in the same Dalitz plot. We parameterise the D 0 → K 0 S π + π − Dalitz plot with the usual variables s + = m 2 (K 0 S π + ) and s − = m 2 (K 0 S π − ), representing the invariant mass-squared of the K 0 S π + and K 0 S π − pair, respectively. If p = (s − , s + ), then p = (s + , s − ), and A To declutter the notation we will omit the superscripts and/or the (p), where there is no risk of ambiguity.
When applying the method presented below in practice it will frequently be necessary to combine results obtained using different conventions for the phase of the CP operator. Switching from our convention with CP |D 0 ⟩ = |D 0 ⟩ to the convention with CP |D 0 ⟩ = −|D 0 ⟩ corresponds to the change δ D → δ D + π.

Measuring
The decay amplitude of a D meson resulting from a B − → DK to a final state f at phase space point p is given by and the corresponding decay rate is (2.4) where N is a normalisation factor. For the CP -conjugate process, with a D from a B + → DK + : For the fits in our validation studies, we follow the widely used practice to re-parameterise the decay rates in terms of the "cartesian" variables This is motivated by the observation that fits in terms of x ± and y ± are statistically better behaved than those in terms of r B , δ B , and γ, which is related to the fact that there is no sensitivity to δ B , γ as r B → 0. The decay rates in terms of the new variables are In a model-dependent approach, A D and A D are derived from a high-statistics amplitude fit to flavour-specific D 0 , D 0 decays. However, the models used to describe the decay amplitude have theoretical shortcomings that make their phase information, which enters via δ D , unreliable. This in turn translates into a systematic uncertainty on γ.
For the model-independent approach described in [4,18], one integrates over regions (bins) of phase space. These bins are defined such that they form CP -conjugate pairs and we label them such that the CP -conjugate of bin i is bin −i. In what follows, we assume the D decays to K 0 S π + π − , although the approach clearly generalises to other decay modes. We define the following parameters related to |A D | 2 , |A D | 2 : In the absence of CP violation in charm, F i = F −i . We also define the following parameters related to the phase-difference between A D and A D : which implies c −i ≡ c i and s −i ≡ −s i . In terms if these quantities, the decay rate B − → DK − with the D decay in phase-space bin i, is given by The CP conjugate process, a B + → DK + decay with the D decay in phase-space bin −i, is Equivalently, these decay rates can be expressed in terms of x ± , y ± : All parameters related to the charm decay, i.e. F i , c i , and s i , can be directly obtained from data, where data from the charm threshold are critical to constraining the parameters related to the phase difference between A D and A D , i.e. c i and s i . It is worth noting, though, that the c i and s i can be obtained alongside γ from (a sufficiently large sample of) B ± → DK ± decays [4]; however, the input from the charm threshold dramatically improves the fit. Measurements of γ in this way, using the D → K 0 S π + π − Dalitz plot we focus on here, have been made by Belle [28,29], Belle & Belle II [30], and LHCb [31][32][33], using input from CLEO-c [34,35] and BES III [36,37]. We will below introduce a new method, also based on exploiting threshold data, that does not require binning, or other forms of integration over phase space as in [22,23]. This is motivated by the aim to maximise the use of information contained in the full two-dimensional Dalitz plot.

Charm threshold data
At the charm threshold, ψ(3770) are produced and decay approximately 50% of the time to a pair of neutral D mesons that we label D 1 , D 2 . The resulting system of D mesons must be C -odd, like the ψ(3770) they originate from. Therefore: (2.16) Let D 1 decay to final state f at phase-space point p and D 2 to final state g at phase-space point q. Then the decay amplitude for this process is We will call g the tag and f the signal (in our feasibility study, f will be K 0 S π + π − ). Depending on the tag, we can distinguish a few important special cases: 1. D 2 → g is a flavour-specific decay such as a semileptonic decay, or a quasi-flavourspecific decay such as D 0 → K − π + (in our feasibility studies, we will ignore the small dilution effect due to the suppressed decay D 0 → K − π + , although this can be taken into account). Then and similarly We refer to these decays as flavour tagged. They provide the same information as flavour-tagged D 0 decays at the B-factories and LHCb, where the flavour is usually identified through the charge of the pion in D * + → D 0 π + and D * − → D 0 π − decays. Measurements of flavour-tagged decays provide |A D (p)| and, in the binned approach, F i . (In recent LHCb analyses, though, F i , and associated efficiency effects, have been obtained from simultaneous fits to B ± → DK ± and B ± → Dπ ± data [33,38].) 2. D 2 → g is a CP -specific decay, either a CP -even decay such as where the expressions for the superpositions of D 0 and D 0 are convention-dependent; we use a convention where CP |D 0 ⟩ = |D 0 ⟩. The corresponding decay amplitudes are 20) and the decay rates: These provide information on cos(δ D ), or, in the binned approach, c i .

(2.24)
This category is therefore sensitive to both sin(δ D ) and cos(δ D ), in contrast to CPtagged decays which are only sensitive to cos(δ D ). In the binned approach, this translates into unique sensitivity to s i .
Measurements of c i and s i for D → K 0 S π + π − have been made at CLEO-c [34,35] and more recently BES III [36,37].
3 Quasi Model Independent unbinned method to measure γ

Basic idea
This section will introduce the quasi model-independent (QMI) method that is the subject of this paper. For concreteness, we will consider the decay D 0 → K 0 S π + π − , although the method generalises to all self-conjugate charm decays, and with some modification also to non-self-conjugate ones.
The statistically most precise way of measuring γ is the model-dependent method. We observe that the magnitude of the amplitude structure across the Dalitz plot of D 0 → K 0 S π + π − is well known from flavour-specific D 0 and D 0 decays analysed by BaBar and Belle [39]. The D 0 , D 0 datasets used in these fits are orders of magnitude larger than the B ± → DK ± datasets available for γ measurements. The fact that the collaborations achieve a decent fit of their models to these datasets implies that we can trust the magnitude of existing amplitude models. However, those models violate unitarity and analyticity, which breaks the connection between magnitude and phase. Consequently, the phases of these models stand on less firm ground. Our approach is therefore to correct the model's phase in a model-independent way, or, more precisely, correct the phase difference between the D 0 amplitude A D (p) and D 0 amplitude A D (p) at each phase space point p. These phase differences are all that matters and in fact it is all we have access to. Measurements of c i and s i by CLEO-c and BES III suggest that the model's phase differences are at least approximately correct [34][35][36][37]. Our approach is therefore to obtain the model-independent phase difference by adding a correcting term δ corr D to the model's phase difference δ model We parameterise δ corr D in a generic way, as a power series in the Dalitz plot parameters. We assume CP conservation in charm decays, such that which also implies This symmetry reduces the number of parameters needed to parameterise δ corr D . We will see below that even if we depart from the assumption that the model's phase differences are approximately correct, such that δ corr D becomes quite sizeable, our method still works. The information that allows us to constrain δ corr D comes, as for the binned methods, predominantly from the charm threshold, although B ± → DK ± decays also contribute.

Constructing the correction to δ D
In order to implement the symmetry relation Eq 3.3, we define the variables Now the symmetry condition from Eq 3.3 becomes δ corr D (z + , z − ) = −δ corr D (z + , −z − ) and can be implemented by only allowing terms with odd powers of z − in the correcting polynomial.
We found that parameterising the phase in terms of Legendre polynomials works well. These are defined for values x ∈ [−1, 1]. We therefore scale z + , z − where the superscripts max and min indicate the maximum and minimum values of the unprimed parameters. The kinematically allowed region of the K 0 S π + π − in terms of the normalised, rotated parameters z ′ + , z ′ − is shown in Fig 1b. The  space where there are no data, with the kinematically-allowed interval in z ′ − varying a lot as a function of z + . We therefore "stretch" the Dalitz plot by replacing z ′ − with The Dalitz plot phase space for the different parametrisations is shown in Fig 1. Alternatives such as a variation of the square Dalitz plot [40] would have achieved a similar outcome. We construct the correcting polynomial, δ corr D , of order N with free parameters C where P n are, in our implementation, n th order Legendre polynomials, although any function with P 2j+1 (z ′′ − ) = −P 2j+1 (−z ′′ − ) would be equally valid. The coefficients C ij are Tag (g) Events CP -even, e.g. Table 1: ψ(3770) → D 1 D 2 , D 1 → K 0 S π + π − , D 2 → g, decays generated for our simulation studies. D 1 , D 2 represent superpositions of D 0 and D 0 , depending on the tag. CP even (odd) tags imply CP -odd (even) D → K 0 S π + π − decays and D 0 (D 0 ) flavour tags imply D 0 (D 0 ) decays to K 0 S π + π − . Sample sizes are taken from [37].
free parameters, determined together with x ± , y ± in a simultaneous fit to ψ(3770) and B ± → DK ± data.

Simulation Studies
We test the algorithm in simulation studies where we generate charm threshold and B ± → DK ± data according to the equations given in Sec 2.2 and 2.3, using the D 0 → K 0 S π + π − model in AmpGen [41]. We ignore backgrounds and efficiency effects in this study.
Using a modified version of AmpGen, we also generate data where we apply a change to the phase of the decay amplitude relative to the model assumed in the fit. We will then show how our new unbinned approach removes the bias this would inflict on the fitted value of γ in a model-dependent approach, and how it does so with improved statistical uncertainty compared to the binned method.

Sample sizes
In our default settings, we simulate charm threshold data according to the event numbers given in BES III's latest D → K 0 S π + π − analysis [37] and B ± → DK ± , D → K 0 S π + π − event yields reported in LHCb's latest measurement of γ in this decay mode [33]. The BES III signal yields for the different tags are shown in Tab 1. For the purpose of this study, we treat D 0 → K − π + and D 0 → K + π − as pure flavour tags and ignore the small contributions from D 0 → K − π + and D 0 → K + π − . LHCb report a total of 12, 533 B ± → DK ± , D → K 0 S π + π − decays. We split this evenly between B − → DK − and B + → DK + , and generate on average 6267 events for each. This implies that, as in LHCb's analysis, our fits are entirely based on the distribution of events within each Dalitz plot, not B + → DK + and B − → DK − event yields integrated across the whole Dalitz plot. We also consider scenarios with larger datasets -by a factor of 100 for B ± → DK ± , D → K 0 S π + π − and by a factor of 10 for charm threshold data.

Amplitudes with modified phases
We generate events based on the nominal amplitude model from BaBar and Belle's joint analysis [39], but with a modified phase difference Note that the expressions for decay rates (2.4, 2.5, 2.7 2.8, 2.21, 2.24) only depend on the phase differences between A D and A D , never the absolute phases of A D or A D themselves. The function f (p) modifies this phase difference. Our approach assures that the magnitudes |A D |, |A D | remain unchanged. We generate events according to three scenarios: 1. no phase modification, f 0 = 0,

double Gaussian modification
The purpose of the error function, erf(x) ≡ 2 √ π x 0 e −t 2 dt, in the definition of f single is to implement the condition f (s + , s − ) = −f (s − , s + ) while providing a smooth transition across the line s + = s − . The parameter values used for scenarios two and three are given in Tab 2. The phase modifications they induce are shown in Fig 2, which shows f single and f double , in radians. It can be seen that scenario two targets the region of the K * resonance, while scenario three has large phase modifications especially in the region of the ρK 0 S − K * π interference. It is worth noting that the phase change relative to the nominal model we consider here is up to ±1 radian, which is not a small shift.

Order by order fits to individual samples
We perform fits with the model-dependent (MD) method and with the quasi model-independent (QMI) method with phase-correction polynomials of order N = 1, 2, . . . , 9, for one sample of each of the phase-modification scenarios. Beyond N = 9, fits converge very slowly due to the large number of parameters. The results are shown in tables 3, 4, and 5, in the format (fit result) − (input value) ± (uncertainty), in units of 10 −2 ; the uncertainty is that estimated by the MINUIT2 [43]-based AmpGen [41] fitter (validated below, in Sec 4.2.2). We see in Tab 3, where there is no phase modification relative to the nominal amplitude model, that the model-dependent and the quasi-model-independent fit both reproduce the input values, within uncertainties. For correction polynomials of order N > 3, the fit results differ slightly between the methods, but far less than the statistical uncertainty. The validation studies below show that this does not lead to a systematic bias. Tables 4 and 5 show that the phase modifications induce a significant bias in x ± , y ± in the model-dependent method, and how the QMI method recovers from it. The changes in the fit results for x ± , y ± between the model-dependent and the higher order QMI fits correspond to changes in the estimated values of r B , δ B , and γ of 0.006, 11 • , and 7 • for the double-Gaussian bias (Tab 5), and of 0.021, 23 • , and 0.3 • for the single-Gaussian bias (Tab 4). An interesting feature is that the uncertainties on x ± and y ± do not appear to be affected significantly by the additional fit parameters. This conclusion is confirmed in the validation studies shown below.
+0.3 ± 0.8 +0.5 ± 1.0 −0.9 ± 0.8 +0.7 ± 1.0 8 +0.3 ± 0.8 +0.5 ± 1.0 −0.9 ± 0.8 +0.7 ± 1.0 9 +0.3 ± 0.8 +0.5 ± 1.0 −0.7 ± 0.8 +0.7 ± 1.0   phase-mod Method  Table 7: Pull results from 100 fits with the quasi-model-independent (QMI) and the modeldependent (MD) method, for each of the three phase-modification scenarios. The QMI fit uses a 6 th order correction polynomial. Results are shown in the format (mean pull) ± (standard deviation). The standard deviation is that of the pull distribution (rather than the uncertainty on the mean). The uncertainty on the mean is 0.1, that on the standard deviation is 0.07. The substantial (and expected) biases observed for the model-dependent method for the fits with phase-modification disappear with the QMI method.

Fits to 100 pseudoexperiments
We generate 100 pseudoexperiments and fit them with the model-dependent method, the binned method, and the QMI method. The QMI method uses a 6 th order correction polynomial. Table 6 shows the mean and standard deviation of the distribution of residuals (i.e. fit result minus truth value) for these fits. Table 7 shows the corresponding value for the pull, which is the residual divided by the uncertainty reported by the fit. The results show that, in the absence of any phase modification, the model-dependent and the QMI method both yield unbiased results with essentially the same uncertainty. The pulls for the QMI method, and for the unbiased data also those for the model-dependent method, show generally good agreement with the expected mean of zero and standard deviation  Table 8: Comparing the QMI method with our implementation of the model-dependent method and the binned method with "optimal" binning (defined in [34]), for the case with no phase modification. The uncertainties given are the average of those reported by the fitter for 100 fits. For the binned fit, we fix c i and s i . In contrast to the QMI results, the uncertainties from the binned fit therefore do not include the effect of the finite sample size at the charm threshold, which leads to an additional uncertainty on γ of 1.2 • [36].
of one. For both methods, the uncertainty the fitter reports on x + seems to be slightly over-estimated. There appears to be a slight under-estimation of the y − uncertainty in the f single configuration. The results with the two phase modification scenarios confirm that the phase modifications induce significant biases in the fit results of the model-dependent method, while the QMI results remain unbiased. Table 8 compares the uncertainties obtained with our new method to those from the model-dependent and the model-independent binned method. Studies in [22] show that the unbinned model-independent method introduced there, which is based on projecting the two-dimensional Dalitz plot onto one dimension, results in a statistical uncertainty on γ between that of the model-dependent and the binned method. The authors of [23] report for their Kolmogorov-Smirnov-inspired unbinned method, for similar simulated event numbers as used in our studies, a statistical uncertainty on γ of ∼ 5 • . However, because of the different values assumed for γ and δ B and differences in the implementation of the amplitude model, comparing the results from [23] with those in Tab 8 is not straightforward.
In our implementation of the binned method, we base the binning on the same amplitude model that we use to generate the simulated data, which should result in a slightly optimistic performance of the binned method. We test all binning schemes defined in [34] and find that the "optimal" binning scheme leads to the best results. In our binned fit, we fix c i and s i to their true value (according to our model), so that the uncertainty on γ for the binned method does not include the contribution from the uncertainty on c i and s i . The sensitivity studies reported in [37] show that, for the "optimal" binning scheme, taking into account the measurement uncertainties on c i and s i leads to an additional uncertainty on γ of 1.2 • . This results in a total uncertainty on γ of 5.1 • ⊕ 1.2 • = 5.2 • , which is, on the same simulated signal data, improved to 4.2 • by the new method introduced, here.

Alternative sample sizes
For the studies above, we used sample sizes corresponding to those reported in recent publications by BES III [36,37] and LHCb [33]. Here we consider possible future datasets that are considerably larger, 10× as large for BES III and 100× as large for LHCb. The results for the QMI method are presented in Tab 10. The results for the model-dependent LHCb σ x + · 10 2 σ y + · 10 2 σ x − · 10 2 σ y − · 10 2 σ γ (  Table 9: Uncertainties on fit parameters for the model-dependent method and the binned method with fixed c i , s i for 1× and 100× the dataset analysed by LHCb in [33]. The uncertainties are the average of those reported by the fitter for fits to 100 pseudoexperiments, generated without backgrounds or detector effects. The statistical uncertainty on the mean of σ γ ranges from 1% to 3% of its value.
Lumi scenario: LHCb BES III σ x + · 10 2 σ y + · 10 2 σ x − · 10 2 σ y − · 10 2 σ γ (  and 100× the dataset analysed by LHCb in [33], and 1× and 10× the dataset analysed by BES III in [36,37]. The uncertainties are the average uncertainty reported by the fitter for ∼ 100 simulated datasets, generated without backgrounds or detector effects. The statistical uncertainty on the mean of σ γ is ∼ 1% of its value. method and the binned method with fixed c i , s i are given in Tab 9. The results for the binned method, with fixed c i , s i , represent the best possible uncertainty on γ that can be reached with this method in the limit of infinitely large threshold data sets, given the 8bin-pair binning scheme and other parameters used in this study. It is possible that a finer binning would improve the uncertainty for the larger data set. The uncertainties for all methods studied scale to a good approximation with 1/ √ N B , where N B is the number of B ± → DK ± , D → K 0 S π + π − events. That this is so for the QMI method is not a priori obvious, given that it depends also on ψ(3770) → DD (i.e. CLEO-c or BES III) data. This suggests that the QMI method is very efficient in extracting information on δ D (s + , s − ) not only from threshold data, but also from B ± → DK ± , D → K 0 S π + π − . For the 1×LHCb scenario, the lack of significant improvement in the uncertainty on γ for the 10× larger BES III sample is consistent with the earlier observation that the QMI method with the 1×LHCb and 1×BES III scenario already achieves effectively the same statistical precision as the model-dependent method (i.e. the best possible for the B ± → DK ± , D → K 0 S π + π − dataset). For the 100×LHCb dataset, the larger BES III sample improves the precision slightly from 0.45 • to 0.43 • , compared to the benchmark of 0.42 • set by the model-dependent method.
Overall, our results indicate that, with the QMI method introduced here, the statistical precision on γ in B ± → DK ± , D → K 0 S π + π − remains close to the optimum defined by the model-dependent method not only with current datasets, but also much larger B ± → DK ± datasets that might become available in the future.

Conclusion
We have introduced a novel unbinned quasi model-idependent (QMI) method for the measurement of γ in B ± → DK ± decays, with input from quantum-correlated charm threshold data. The method uses a polynomial to correct the phase of the D meson's decay amplitude model in an unbinned, model-independent way. We studied the performance of the method using simulated B ± → DK ± , D → K 0 S π + π − and ψ(3770) → DD signal events. The method produces unbiased results for cases where discrepancies between the assumed amplitude model and the true model produce large biases in a model-dependent measurement. For realistic current and plausible future sample sizes, the method achieves a statistical precision on γ that is effectively the same as the optimum defined by the model-dependent method, without suffering from the systematic uncertainty associated with the amplitude model. The statistical uncertainty is significantly better than that of the binned model-independent method currently in use.
We expect that the QMI method will also be beneficial in the study of charm mixing and the study of phases of decay amplitudes across the Dalitz plot.