The Belle Collaboration

In this letter we report the first simultaneous determination of the absolute value of the Cabibbo-Kobayashi-Maskawa (CKM) matrix element $V_{ub}$ using inclusive and exclusive decays. Exclusive determinations of $|V_{ub}|$ focus on reconstructing explicit final states such as $B\to\pi\ell\bar{\nu}_{\ell}$ Amhis et al. (2023), $\Lambda_{b}\to p\mu\bar{\nu}_{\mu}$ Aaij et al. (2015), or $B_{s}\to K\,\mu\bar{\nu}_{\mu}$ Aaij et al. (2021), whereas inclusive determinations study $B$ meson decays undergoing $b\to u\ell\bar{\nu}_{\ell}$ transitions and are indiscriminate of the $u\to X_{u}$ hadronization process. The world averages of either method are only marginally compatible Amhis et al. (2023),

with a ratio of $|V_{ub}^{\mathrm{excl.}}|/|V_{ub}^{\mathrm{incl.}}|=0.84\pm 0.04$, which deviates 3.7 standard deviations from unity. The underlying reason for this tension is unknown. New physics explanations are challenging (see e.g. Refs. Crivellin (2010); Enomoto and Tanaka (2015); Bernlochner et al. (2014); Umeeda (2022)), leading to some to speculate the existence of until now unaccounted systematic effects Zyla et al. (2020a). This motivates the simultaneous determination in a single analysis, in which $B\to\pi\ell\bar{\nu}_{\ell}$ and the $B\to X_{u}\,\ell\,\bar{\nu}_{\ell}$ rates can be simultaneously extracted and systematic effects can be correlated.

The presented measurement of inclusive and exclusive $b\to u\ell\bar{\nu}_{\ell}$ decays uses the same collision events and a similar analysis strategy as Refs. Cao et al. (2021a, b). Charmless semileptonic decays are reconstructed by relying on the complete reconstruction of the second $B$ meson in the $e^{+}e^{-}\to\Upsilon(4S)\to B\bar{B}$ process. This approach allows for the direct reconstruction of the $X_{u}$ system of the $B\to X_{u}\,\ell\,\bar{\nu}_{\ell}$ process. Specifically, the four-momentum transfer squared, $q^{2}=\left(p_{B}-p_{X_{u}}\right)^{2}$, and the number of charged pion candidates of the $X_{u}$ system, $N_{\pi^{\pm}}$, can be reconstructed. This allows for the separation of $B^{+}\to\pi^{0}\,\ell^{+}\,\nu_{\ell}$ and $B^{0}\to\pi^{-}\,\ell^{+}\,\nu_{\ell}$ from other $B\to X_{u}\,\ell\,\bar{\nu}_{\ell}$ decays. The main background in the measurement stems from the much more abundant $B\to X_{c}\,\ell\,\bar{\nu}_{\ell}$ decays and a multivariate suppression strategy is used to reduce this and other background processes. Charge conjugation is implied throughout. The inclusive $B\to X_{u}\,\ell\,\bar{\nu}_{\ell}$ branching fraction is defined as the average branching fraction of $B^{+}$ and $B^{0}$ meson decays. Furthermore, we denote $\ell=e,\mu$, and use natural units: $\hbar=c=1$.

We analyze $(772\pm 10)\times 10^{6}$ $B$ meson pairs recorded at the $\Upsilon(4S)$ resonance energy and $79\text{\,}\mathrm{f}\mathrm{b}^{-1}$ of collision events recorded $60\text{\,}\mathrm{M}\mathrm{e}\mathrm{V}$ below the $\Upsilon(4S)$ peak. Both data sets were recorded at the KEKB $e^{+}e^{-}$ collider Kurokawa and Kikutani (2003) by the Belle detector. Belle is a large-solid-angle magnetic spectrometer. A detailed description of its performance and subdetectors can be found in Ref. (Abashian et al., 2002). The particle identification and selection criteria are the same as in Ref. Cao et al. (2021a).

Monte Carlo (MC) samples of $B$ meson decays and continuum processes ($e^{+}e^{-}\to q\bar{q}$ with $q=u,d,s,c$) are simulated using the EvtGen generator (Lange, 2001). The normalization of continuum events is calibrated with the measured off-resonance data. A detailed description of all samples and decay models is given in Ref. Cao et al. (2021a). The simulated samples are used for background subtraction and to correct for detector resolution, selection, and acceptance effects. The used sample sizes correspond to approximately ten and five times, respectively, the Belle collision data for the $B$ meson production and continuum processes.

Semileptonic $B\to X_{u}\,\ell\,\bar{\nu}_{\ell}$ decays are simulated as a mixture of specific exclusive modes and nonresonant contributions using a “hybrid” approach Ramirez et al. (1990); Prim et al. (2020); Prim (2020): the triple differential rate of inclusive and exclusive predictions are combined such that the partial rates of the inclusive prediction are recovered. This is achieved by assigning weights to the inclusive contribution as a function of the generator-level $q^{2}$, $E_{\ell}^{B}$, and $M_{X}$. Here $E_{\ell}^{B}$ and $M_{X}$ denote the energy of the lepton in the signal $B$ rest frame and the invariant mass of the $X_{u}$ system produced in the $B\to X_{u}\,\ell\,\bar{\nu}_{\ell}$ decay. For the inclusive contribution, we use two different calculations: the De Fazio and Neubert (DFN) model De Fazio and Neubert (1999) (with $m_{b}^{\text{KN}}=(4.66\pm 0.04)\,\mathrm{GeV}$, $a^{\text{KN}}=1.3\pm 0.5$) and the Bosch-Lange-Neubert-Paz (BLNP) model Lange et al. (2005) (with $m_{b}^{\mathrm{SF}}=4.61\,\mathrm{GeV}$, $\mu_{\pi}^{2\,\text{SF}}=0.20\,\mathrm{GeV}^{2}$). The difference between the two models is treated as a systematic uncertainty. The simulated inclusive $B\to X_{u}\,\ell\,\bar{\nu}_{\ell}$ events are hadronized with the JETSET algorithm T. Sjöstrand (1994) into final states with two or more mesons. We study two different tunes of the fragmentation parameters and assign their difference as a systematic uncertainty. The exclusive contributions are modeled as follows: $B\to\pi\,\ell\,\bar{\nu}_{\ell}$ decays are modeled using the Bourrely-Caprini-Lellouch (BCL) form factor parameterization (Bourrely et al., 2009); $B\to\rho\,\ell\,\bar{\nu}_{\ell}$ and $B\to\omega\,\ell\,\bar{\nu}_{\ell}$ decays are modeled using the Bharucha-Straub-Zwicky (BSZ) form factors (Bharucha et al., 2016) from the fit of Ref. (Bernlochner et al., 2021) to light-cone sum rule (LCSR) predictions (Bharucha et al., 2016) and the measurements of Refs. Sibidanov et al. (2013a); Lees et al. (2013); del Amo Sanchez et al. (2011a); $B\to\eta\,\ell\,\bar{\nu}_{\ell}$ and $B\to\eta^{\prime}\,\ell\,\bar{\nu}_{\ell}$ are modeled using pole form factors obtained from fits to LCSR (Duplancic and Melic, 2015). For the branching fractions the world averages from Ref. (Zyla et al., 2020b) are used.

Semileptonic $B\to X_{c}\,\ell\,\bar{\nu}_{\ell}$ decays are dominated by $B\to D\,\ell\,\bar{\nu}_{\ell}$ and $B\to D^{*}\,\ell\,\bar{\nu}_{\ell}$ decays. We simulate them with the form factors of Refs. Boyd et al. (1995); Grinstein and Kobach (2017); Bigi et al. (2017) and values determined by the measurements of Refs. Glattauer et al. (2016); Waheed et al. (2019). Other $B\to X_{c}\,\ell\,\bar{\nu}_{\ell}$ decays are simulated as a mixture of resonant and nonresonant modes, using the parameterization of Ref. Bernlochner and Ligeti (2017) for the modeling of $B\to D^{**}\,\ell\,\bar{\nu}_{\ell}$ form factors. The known difference between inclusive and the sum of measured exclusive $B\to X_{c}\,\ell\,\bar{\nu}_{\ell}$ decays is simulated with $B\to D^{(*)}\,\eta\,\ell^{+}\nu_{\ell}$ decays.

We reconstruct $e^{+}e^{-}$ collision events with the multivariate tagging algorithm of Ref. Feindt et al. (2011). The algorithm uses a hierarchical approach utilizing neural networks to fully reconstruct one of the two $B$ mesons in hadronic final states (labeled as $B_{\mathrm{tag}}$). The $B_{\mathrm{tag}}$ reconstruction efficiency is calibrated using $B\to X_{c}\,\ell\,\bar{\nu}_{\ell}$ decays following the prescription outlined in Cao et al. (2021a). The identified final state particles forming the $B_{\mathrm{tag}}$ are masked and $b\to u\ell\bar{\nu}_{\ell}$ signal candidates are reconstructed by identifying an electron or muon candidate in the events, requiring $E_{\ell}^{B}=|\mathbb{p}_{\ell}^{B}|>1\,\mathrm{GeV}$ as measured in the signal $B$ rest frame. To reject background from the much more abundant $B\to X_{c}\,\ell\,\bar{\nu}_{\ell}$ decays, eleven distinguishing features are combined into a single discriminant using boosted decision trees (BDTs) and utilizing the implementation of Ref. Chen and Guestrin (2016). The most discriminating training features are the reconstructed neutrino mass, $M_{\mathrm{miss}}^{2}$, the vertex fit probability of the decay vertex between the hadronic system $X$ and the signal lepton $\ell$, and the number of identified $K^{\pm}$ and $K_{S}^{0}$ in the $X$ system. Same as in Cao et al. (2021a), we select a working point that corresponds to a signal efficiency of 18.5%, which rejects 98.7% of $B\to X_{c}\,\ell\,\bar{\nu}_{\ell}$ decays, defined with respect to all events after the $B_{\mathrm{tag}}$ selection. To test the modeling of $B\to X_{c}\,\ell\,\bar{\nu}_{\ell}$ and other backgrounds in the extraction variables, $q^{2}$ and $N_{\pi^{\pm}}$, we also utilize the events failing the BDT selection and find good agreement sup . We further separate events by the reconstructed $M_{X}$, categorizing $M_{X}<1.7\,\mathrm{GeV}$ into five $q^{2}$ bins ranging in $[0,26.4]\,\mathrm{GeV}^{2}$ as a function of the $N_{\pi^{\pm}}$ multiplicity for the interval of $[0,1,2,\geq 3]$. Events with $M_{X}\geq 1.7\,\mathrm{GeV}$ are analyzed only in bins of $N_{\pi^{\pm}}$ as they are dominated by background. To enhance the $B\to\pi\,\ell\,\bar{\nu}_{\ell}$ purity in the low-$M_{X}$ $N_{\pi^{\pm}}=0$ and $N_{\pi^{\pm}}=1$ events, we apply a selection on the thrust of 0.92 and 0.85, respectively. It is defined by $\max_{|\mathbb{n}|=1}\left(\sum_{i}|\mathbb{p_{i}}\cdot\mathbb{n}|/\sum_{i}|\mathbb{p_{i}}|\right)$, when summing over the neutral and charged constituents of the reconstructed $X$ system in the center of mass frame. For $B\to\pi\,\ell\,\bar{\nu}_{\ell}$ events, we expect a more collimated $X_{u}$ system than for $B\to X_{c}\,\ell\,\bar{\nu}_{\ell}$ and other $B\to X_{u}\,\ell\,\bar{\nu}_{\ell}$ processes, resulting in a higher thrust value.

The $q^{2}:N_{\pi^{\pm}}$ bins and the $M_{X}\geq 1.7\,\mathrm{GeV}$ $N_{\pi^{\pm}}$ distribution are analyzed using a simultaneous likelihood fit, which incorporates floating parameters for the modeling of the $B\to\pi\,\ell\,\bar{\nu}_{\ell}$ form factor, the binned templates, and systematic uncertainties as nuisance parameters. Specifically, the shape of $B\to\pi\,\ell\,\bar{\nu}_{\ell}$ template is linked to the form factors by correcting the efficiency and acceptance effects. The fit components we probe are the normalizations of $B\to\pi\,\ell\,\bar{\nu}_{\ell}$ decays, other $B\to X_{u}\,\ell\,\bar{\nu}_{\ell}$ signal decays, and of background events dominated by $B\to X_{c}\,\ell\,\bar{\nu}_{\ell}$ decays. The $f_{+}$ and $f_{0}$ form factors describing the $B\to\pi\,\ell\,\bar{\nu}_{\ell}$ decay dynamics are parameterized with expansion coefficients $a_{n}^{+}$ and $a_{n}^{0}$ using the BCL expansion,

at expansion order $N^{+}=N^{0}=3$ in the conformal variable $z=z(q^{2})$ (Bourrely et al., 2009; Aoki et al., 2022), and $a^{0}_{2}$ is expressed by the remaining coefficients to keep the kinematical constraint $f_{+}(0)=f_{0}(0)$. We constrain the expansion coefficients to the lattice QCD (LQCD) values of Ref. Aoki et al. (2022), combining LQCD calculations from several groups Bailey et al. (2015); Flynn et al. (2015). Note that the measured distributions have no sensitivity for $f_{0}$ and we thus neglect its effects in the decay rate. The inclusion of the $f_{0}$ expansion coefficients, however, reduces uncertainties on the $B\to\pi\,\ell\,\bar{\nu}_{\ell}$ rate through the correlation to the $f_{+}$ shape. In order to utilize the full experimental knowledge of the $B\to\pi$ form factors to date, we constrain its shape to the combined lattice QCD and experimental information of Refs. del Amo Sanchez et al. (2011b); Lees et al. (2012); Ha et al. (2011); Sibidanov et al. (2013b). The fit scenario with only lattice QCD constraints is studied for a standalone comparison with other experimental results.

We consider additive and multiplicative systematic uncertainties in the likelihood fit by adding bin-wise nuisance parameters for each template. The parameters are constrained to a multinormal Gaussian distribution with a covariance reflecting the sum of all considered systematic effects, and the correlation structure between templates from common sources is taken into account. This includes detector and reconstruction related uncertainties, such as the tracking efficiency for low and high momentum tracks, particle identification efficiency uncertainties, and the calibration of the $B_{\mathrm{tag}}$ reconstruction efficiency. We further consider uncertainties on the $B\to X_{u}\,\ell\,\bar{\nu}_{\ell}$ and $B\to X_{c}\,\ell\,\bar{\nu}_{\ell}$ shapes from the form factors, non-perturbative parameters, and their compositions. The $u\to X_{u}$ fragmentation uncertainties are evaluated by changing the default Belle tune of fragmentation parameters to the values used in Ref. Aaij et al. (2014). We further vary the $s\bar{s}$-production rate $\gamma_{s}=0.30\pm 0.09$, spanning the range of Refs. Althoff et al. (1985); Bartel et al. (1983). The largest uncertainties on the exclusive branching fraction measurements are from the calibration of the tagging efficiency ($\pm 4.1\%$) and the $B\to X_{u}\,\ell\,\bar{\nu}_{\ell}$ modeling ($\pm 3.5\%$). The largest uncertainties on the inclusive branching fraction measurement are from the $B\to X_{u}\,\ell\,\bar{\nu}_{\ell}$ ($\pm 10.9\%$) modeling and the $u\to X_{u}$ fragmentation ($\pm 5.3\%$). The uncertainties of the modeling of the $B\to X_{c}\,\ell\,\bar{\nu}_{\ell}$ background are $\pm 1.2\%$ and $\pm 2.8\%$ for the $B\to\pi\,\ell\,\bar{\nu}_{\ell}$ and $B\to X_{u}\,\ell\,\bar{\nu}_{\ell}$ branching fractions, respectively.

Figure 1 shows the $q^{2}:N_{\pi^{\pm}}$ distribution of the signal region after the fit and with only using LQCD information: $B^{+}\to\pi^{0}\,\ell^{+}\,\nu_{\ell}$ and $B^{0}\to\pi^{-}\,\ell^{+}\,\nu_{\ell}$ events are aggregated in the $N_{\pi^{+}}=0$ and $N_{\pi^{+}}=1$ bins, respectively, whereas contributions from other $B\to X_{u}\,\ell\,\bar{\nu}_{\ell}$ processes are in all multiplicity bins. The high $M_{X}$ bins constrain the $B\to X_{c}\,\ell\,\bar{\nu}_{\ell}$ and other background contributions. We use the isospin relation and $B^{0}/B^{+}$ lifetime ratio to link the yields of $B^{+}\to\pi^{0}\,\ell^{+}\,\nu_{\ell}$ and $B^{0}\to\pi^{-}\,\ell^{+}\,\nu_{\ell}$. The fit has a $\chi^{2}$ of $13.8$ with $21$ degrees of freedom, corresponding to a p-value of $88\%$. The measured $B^{+}\to\pi^{0}\,\ell^{+}\,\nu_{\ell}$ and $B^{0}\to\pi^{-}\,\ell^{+}\,\nu_{\ell}$ yields are corrected for efficiency effects to determine the corresponding branching fractions $\mathcal{B}$. The measured inclusive yield is calculated from the sum of $B^{+}\to\pi^{0}\,\ell^{+}\,\nu_{\ell}$, $B^{0}\to\pi^{-}\,\ell^{+}\,\nu_{\ell}$, and other $B\to X_{u}\,\ell\,\bar{\nu}_{\ell}$ events and unfolded to correspond to a partial branching fraction $\Delta\mathcal{B}$ with $E_{\ell}^{B}>1.0\,\mathrm{GeV}$, also correcting for the effect of final state radiation photons. We find

with the errors denoting statistical and systematic uncertainties and we used the isospin relation between $B^{-}\to\pi^{0}\ell^{-}\bar{\nu}_{\ell}$ and $\overline{B}^{0}\to\pi^{+}\ell^{-}\bar{\nu}_{\ell}$ to link both branching fractions. The recovered branching fraction for $\overline{B}^{0}\to\pi^{+}\ell^{-}\bar{\nu}_{\ell}$ is compatible with the world average of $\mathcal{B}(\overline{B}^{0}\to\pi^{+}\ell^{-}\bar{\nu}_{\ell})=\left(1.50\pm 0.06\right)\times 10^{-4}$  Amhis et al. (2023). The correlation between the exclusive and inclusive branching fractions is $\rho=0.11$. Using calculations for the inclusive partial rate and the fitted form factor parameters, we can determine values for $|V_{ub}|$. As our baseline we use the GGOU Gambino et al. (2007) calculation for the inclusive partial rate with $E_{\ell}^{B}>1.0\,\mathrm{GeV}$ ($\Delta\Gamma/|V_{ub}|^{2}=58.5\pm 2.7$ ps^-1), but other calculations result in similar values for inclusive $|V_{ub}|$. We find

for exclusive and inclusive $|V_{ub}|$ with the uncertainties denoting the statistical error, systematic error, and error from theory (either from LQCD or the inclusive calculation). The correlation between the exclusive and inclusive $|V_{ub}|$ is $\rho=0.07$. The determined value for inclusive $|V_{ub}|$ is compatible with the determination of Ref. Cao et al. (2021a). For the ratio of inclusive and exclusive $V_{ub}$ values, we find

which is compatible with the SM expectation of unity. The value is higher and compatible with the current world average of $|V_{ub}^{\mathrm{excl.}}|/|V_{ub}^{\mathrm{incl.}}|=0.84\pm 0.04$ Amhis et al. (2023) within 1.5 standard deviations. Fig. 2 (top) compares the measured individual values with the SM expectation and the current world average. We also test what happens if we relax the isospin relation between $B^{-}\to\pi^{0}\ell^{-}\bar{\nu}_{\ell}$ (red ellipse) and $\overline{B}^{0}\to\pi^{+}\ell^{-}\bar{\nu}_{\ell}$ (blue) branching fractions and find compatible results for exclusive and inclusive $|V_{ub}|$, as well as for the exclusive $|V_{ub}|$ values.

In the nominal result, we utilize the full theoretical and experimental knowledge of the $B\to\pi\,\ell\,\bar{\nu}_{\ell}$ form factor, combining shape information from the measured $q^{2}$ spectrum with LQCD predictions, as provided by Ref. Aoki et al. (2022). The determined (partial) branching fractions in this scenario are

with a correlation of $\rho=0.12$ between inclusive and exclusive branching fractions and assuming isospin relation. This fit leads to a more precise value of $|V_{ub}|$ from $B\to\pi\,\ell\,\bar{\nu}_{\ell}$ and we find with the same inclusive calculation

with a correlation $\rho=0.11$ and a ratio of

compatible with the world average within 1.2 standard deviations. Fig. 2 (bottom) compares the obtained values and we also find good agreement between the isospin conjugate exclusive values of $|V_{ub}|$. Figure 3 compares the fitted $q^{2}$ spectra of the differential rate of $\overline{B}^{0}\to\pi^{+}\ell^{-}\bar{\nu}_{\ell}$ for both fit scenarios as well as for the LQCD input Aoki et al. (2022). The inclusion of the full experimental and theoretical knowledge leads to a higher rate at low $q^{2}$.

In summary, we presented the first simultaneous determination of inclusive and exclusive $|V_{ub}|$ within a single analysis. In the ratio of both $|V_{ub}|$ values many systematic uncertainties such as the tagging calibration or the lepton identification uncertainties cancel and one can directly test the SM expectation of unity. We recover ratios that are compatible with this expectation, but 1.5 standard deviations higher than the ratio of the current world averages of inclusive and exclusive $|V_{ub}|$. This tension is reduced to 1.2 standard deviations when including the constraint based on the full theoretical and experimental knowledge of the $B\to\pi\,\ell\,\bar{\nu}_{\ell}$ form factor shape. We average our inclusive and exclusive values from both approaches using LQCD or LQCD and additional experimental information and find,

respectively. These values can be compared with the expectation from CKM unitarity of Ref. Charles et al. (2005) of $|V_{ub}^{\mathrm{CKM}}|=(3.64\pm 0.07)\times 10^{-3}$ and are compatible within 1.2 and 0.8 standard deviations, respectively. The applied approach of simultaneously fitting $q^{2}$ and the number of charged pions in the $X_{u}$ system will benefit from the large anticipated data set of Belle II. Additional fit scenarios and inclusive $|V_{ub}|$ values from other theory calculations of the partial rate are provided in the supplemental material sup .

This work, based on data collected using the Belle detector, which was operated until June 2010, was supported by the Ministry of Education, Culture, Sports, Science, and Technology (MEXT) of Japan, the Japan Society for the Promotion of Science (JSPS), and the Tau-Lepton Physics Research Center of Nagoya University; the German Research Foundation (DFG) Emmy-Noether Grant No. BE 6075/1-1; the Helmholtz W2/W3-116 grant; the Australian Research Council including grants DP180102629, DP170102389, DP170102204, DE220100462, DP150103061, FT130100303; Austrian Federal Ministry of Education, Science and Research (FWF) and FWF Austrian Science Fund No. P 31361-N36; the National Natural Science Foundation of China under Contracts No. 11675166, No. 11705209; No. 11975076; No. 12135005; No. 12175041; No. 12161141008; Key Research Program of Frontier Sciences, Chinese Academy of Sciences (CAS), Grant No. QYZDJ-SSW-SLH011; Project ZR2022JQ02 supported by Shandong Provincial Natural Science Foundation; the Ministry of Education, Youth and Sports of the Czech Republic under Contract No. LTT17020; the Czech Science Foundation Grant No. 22-18469S; Horizon 2020 ERC Advanced Grant No. 884719 and ERC Starting Grant No. 947006 “InterLeptons” (European Union); the Carl Zeiss Foundation, the Deutsche Forschungsgemeinschaft, the Excellence Cluster Universe, and the VolkswagenStiftung; the Department of Atomic Energy (Project Identification No. RTI 4002) and the Department of Science and Technology of India; the Istituto Nazionale di Fisica Nucleare of Italy; National Research Foundation (NRF) of Korea Grant Nos. 2016R1D1A1B02012900, 2018R1A2B3003643, 2018R1A6A1A06024970, RS202200197659, 2019R1I1A3A01058933, 2021R1A6A1A03043957, 2021R1F1A1060423, 2021R1F1A1064008, 2022R1A2C1003993; Radiation Science Research Institute, Foreign Large-size Research Facility Application Supporting project, the Global Science Experimental Data Hub Center of the Korea Institute of Science and Technology Information and KREONET/GLORIAD; the Polish Ministry of Science and Higher Education and the National Science Center; the Ministry of Science and Higher Education of the Russian Federation, Agreement 14.W03.31.0026, and the HSE University Basic Research Program, Moscow; University of Tabuk research grants S-1440-0321, S-0256-1438, and S-0280-1439 (Saudi Arabia); the Slovenian Research Agency Grant Nos. J1-9124 and P1-0135; Ikerbasque, Basque Foundation for Science, Spain; the Swiss National Science Foundation; the Ministry of Education and the Ministry of Science and Technology of Taiwan; and the United States Department of Energy and the National Science Foundation. These acknowledgements are not to be interpreted as an endorsement of any statement made by any of our institutes, funding agencies, governments, or their representatives. We thank the KEKB group for the excellent operation of the accelerator; the KEK cryogenics group for the efficient operation of the solenoid; and the KEK computer group and the Pacific Northwest National Laboratory (PNNL) Environmental Molecular Sciences Laboratory (EMSL) computing group for strong computing support; and the National Institute of Informatics, and Science Information NETwork 6 (SINET6) for valuable network support. We are indebted to Alexander Ermakov for his pioneering work on the subject. We thank Frank Tackmann, Zoltan Ligeti, and Dean Robinson for discussions about the content of this manuscript.