Test of lepton flavor universality with measurements of $R(D^{+})$ and $R(D^{*+})$ using semileptonic $B$ tagging at the Belle II experiment

A fundamental property of the Standard Model (SM) of particle physics is the universality of the electroweak gauge couplings to the three fermion generations. In the lepton sector, this universality results in an accidental symmetry of the lepton flavors that is only broken by the Higgs-Yukawa interaction. One key consequence is that physical processes involving charged leptons feature lepton flavor universality (LFU), an approximate symmetry of lepton flavor among physical observables, only broken by charged lepton mass terms emerging from the non-zero vacuum expectation value of the Higgs field. An observation of lepton flavor universality violation would therefore be a clear signature of physics beyond the SM [1].

In this paper, we test LFU using semitauonic $b\to c\tau\bar{\nu}_{\tau}$ decays by measuring the ratios,¹¹1Charge conjugation is implied throughout this paper.

with $\ell=e,\mu$ and $D^{(*)}$ denoting either $D^{(*)+}$ or $D^{(*)0}$. The first equality follows from the assumption of isospin symmetry. Predictions of these ratios are independent of the magnitude of the Cabibbo-Kobayashi-Maskawa (CKM) matrix element $V_{cb}$ and, to some extent, of the parameterization of hadronic matrix elements, reaching a precision of 1–2%  [2, 3, 4, 5, 6, 7, 8, 9, 10]. Experimentally, measurements of the ratios in Eq. 1 are preferred to measurements of absolute branching fractions, as efficiency-related systematic uncertainties largely cancel. Furthermore, the simultaneous measurement of both $\mathcal{R}(D)$ and $\mathcal{R}(D^{*})$ is useful, as $\overline{B}{}\to D^{*}\tau^{-}\bar{\nu}_{\tau}$ decays are an important background in the reconstruction of $\overline{B}{}\to D\tau^{-}\bar{\nu}_{\tau}$ decays.

Several experiments previously reported measurements of these or similar ratios [11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22]. Belle II reported $\mathcal{R}(D^{*})=0.262^{+0.041}_{-0.039}(\mathrm{stat})^{+0.035}_{-0.032}(\mathrm{syst})$ [23] and the inclusive ratio $\mathcal{R}(X_{\tau/\ell})=0.228\pm 0.016\mathrm{(stat)}\pm 0.036\mathrm{(syst)}$ [24]. A combination of these measurements is reported in Ref. [25] and achieves a precision of 8% for $\mathcal{R}(D)$ and 4% for $\mathcal{R}(D^{*})$, with values of $\mathcal{R}(D)=0.342\pm 0.026$ and $\mathcal{R}(D^{*})=0.286\pm 0.012$. Both exceed the SM expectations of $\mathcal{R}(D)=0.296\pm 0.004$ and $\mathcal{R}(D^{*})=0.254\pm 0.005$ [2, 3, 4, 5, 6, 7, 8, 9, 10] with a significance of 3.1 standard deviations.

In this analysis, we study $B^{0}\overline{B}{}^{0}$ pairs produced in $\Upsilon(4S)$ decays. One $B^{0}$ or $\overline{B}{}^{0}$ meson is reconstructed in a semileptonic decay using the hierarchical reconstruction algorithm from Ref. [26], referred to as $B_{\mathrm{tag}}$. This $B_{\mathrm{tag}}$ is then combined with an oppositely flavored semitauonic or semileptonic decay candidate, which we define as $B_{\mathrm{sig}}$. The selected events are independent of those analyzed in Refs. [24, 23], as the latter relied on the reconstruction of hadronically decaying $B_{\mathrm{tag}}$ mesons. The analysis of $B^{+}B^{-}$ pairs is deferred to future work, as the reconstruction of the isospin-conjugate $B^{-}$ signal decays necessitates the efficient identification of $D^{*0}\to D^{0}\pi^{0}$ decays, which suffer from an increased combinatorial background.

We reconstruct signal candidates from the $D^{+}\ell^{-}$ and $D^{*+}\ell^{-}$ final states (with $\ell=e,\,\mu$), which can originate from either semitauonic decays ($\overline{B}{}^{0}\to D^{(*)+}\tau^{-}\bar{\nu}_{\tau}$ with $\tau^{-}\to\ell^{-}\overline{\nu}_{\ell}\nu_{\tau}$) or semileptonic decays ($\overline{B}{}^{0}\to D^{(*)+}\ell^{-}\bar{\nu}_{\ell}$). These processes differ in the number of neutrinos and are thus distinguishable through kinematic properties. Since both processes produce the same visible final state in the detector, several experimental systematic uncertainties cancel in the measurement of their ratio.

We distinguish $B_{\mathrm{sig}}$ candidates originating from semitauonic, semileptonic, and background sources using multivariate classifiers trained on the kinematic properties of both the $B_{\mathrm{sig}}$ and $B_{\mathrm{tag}}$ candidates. A binned maximum likelihood fit is then performed to measure the relative contribution from each source, allowing a direct determination of $\mathcal{R}(D^{(*)})$.

The remainder of this paper is structured as follows: Sections II and III provide an overview of the Belle II detector, the analyzed data set, and the simulated samples. Section IV summarizes the tag and signal reconstruction, while Section V describes the employed multivariate selection. Section VI details the fitting procedure, and Section VII discusses the systematic uncertainties affecting the measurement. Section VIII presents our findings and consisitency checks, and Section IX provides our conclusions.

The analysis uses Belle II data collected at SuperKEKB [27] from 2019 to 2022 at a center-of-mass energy of 10.58 GeV,²²2Natural units ($c=\hbar=1$) are used throughout this paper. corresponding to the $\Upsilon(4S)$ resonance, having an integrated luminosity of 365 $\text{fb}^{-1}$. The sample contains an estimated $(387\pm 6)\times 10^{6}$ $B\overline{B}{}$ events. Additionally, 42.3 fb^-1 of off-resonance data at 10.52 GeV is used to study $q\bar{q}$ ($q=u,d,s,c$) background.

The Belle II detector [28] is an upgraded version of Belle [29] with enhanced particle reconstruction and identification. Its subdetectors are arranged cylindrically around the interaction point (IP), which is enclosed by a $1\text{\,}\mathrm{cm}$ beryllium beam pipe. The pixel detector (PXD) consists of two layers, with the first fully instrumented and the second partially completed. The PXD is surrounded by a four-layer double-sided silicon-strip detector (SVD), and both detectors are used to reconstruct decay vertices with high precision. Surrounding these detectors is the central drift chamber (CDC), which provides three-dimensional tracking and specific ionization ($\mathrm{d}{E}/\mathrm{d}{x}$) measurements.

Outside the CDC, the time-of-propagation (TOP) and aerogel ring-imaging Cherenkov (ARICH) detectors provide particle identification in the barrel and forward endcap regions, respectively.

The electromagnetic calorimeter (ECL), consisting of a $3\text{\,}\mathrm{m}$ barrel and annular endcaps, is located outside the TOP and within a $1.5\text{\,}\mathrm{T}$ superconducting solenoid. The $K^{0}_{L}$ and muon detector (KLM), situated outside the solenoid, is composed of iron plates interleaved with active detector elements.

Particle candidates are constructed and identified using the information from various detector systems.

Charged particle candidates (tracks) are reconstructed by the vertex and tracking systems, and identified based on information from the outer detectors. In particular, muons with sufficiently high momentum will traverse the KLM, while other charged particles are absorbed. In contrast, electrons deposit nearly all of their energy in the ECL.

Photon candidates consist of ECL clusters that are not consistent with extrapolations of charged tracks. Minimum energy selections are necessary to reject clusters from beam-induced background photons.

Monte Carlo (MC) samples are used to determine reconstruction efficiencies and acceptance effects as well as to estimate background contamination and to train multivariate classifiers. The $B$ decays are simulated using the EvtGen generator [30]. The simulation of $e^{+}e^{-}\to q\bar{q}$ continuum processes is carried out with KKMC [31] and PYTHIA8 [32]. Electromagnetic final-state radiation is simulated using PHOTOS [33] for all charged final-state particles. Interactions of particles with the detector are simulated using GEANT4 [34]. The simulated samples contain the equivalent of 2.8 ab^-1 of $B\overline{B}{}$ and continuum processes. The $B\overline{B}{}$ events are simulated with equal fractions of neutral and charged $B$ mesons. An additional sample of $176\times 10^{6}$ $\overline{B}{}^{0}\to D^{(*)+}\tau^{-}\bar{\nu}_{\tau}$ decays with $\tau^{-}\to\ell^{-}\overline{\nu}_{\ell}\nu_{\tau}$ is used, corresponding to an effective sample size of $8.5$ ab^-1.

The signal decays $\overline{B}{}^{0}\to D^{(*)+}\tau^{-}\bar{\nu_{\tau}}$ and $\overline{B}{}^{0}\to D^{(*)+}\ell^{-}\bar{\nu_{\ell}}$ are modeled using the form factors from Ref. [6], with parameter values obtained from a fit to the measurements in Refs. [35, 36]. To incorporate this form factor model into the Monte Carlo simulation, the HAMMER software package [37] is used to compute and apply event-by-event weights. For the branching fractions isospin-averaged values of Ref. [38] are used.

The decays $\overline{B}{}^{0}\to D^{**}\tau^{-}\overline{\nu}_{\tau}$ and $\overline{B}{}^{0}\to D^{**}\ell\overline{\nu}_{\ell}$, where $D^{**}=\{D_{0}^{*+},D_{1}^{\prime+},D_{1}^{+},D_{2}^{*+}\}$, are modeled using the heavy-quark-symmetry-based form factors proposed in Refs. [39, 40], with $D^{**}$ masses and widths taken from Ref. [41]. For the $\overline{B}{}\to D^{**}\ell\overline{\nu}_{\ell}$ branching fractions, we adopt the values from Ref. [38] to account for missing isospin-conjugated and other established decay modes observed in studies of $B$ decays into fully hadronic final states, following the approach outlined in Ref. [39].

The difference between the inclusive semileptonic branching fraction and the sum of exclusive semileptonic $B$ decays (the so-called “gap”) is accounted for by using a dedicated sample of $\overline{B}{}\to D^{(*)}(\pi\pi/\eta)\ell\overline{\nu}_{\ell}$ decays. These are simulated using broad $D^{**}$ contributions that do not form distinct resonance peaks in their invariant mass distributions, resulting in a smooth continuum across phase space. The heavy-quark-symmetry-based form factors of Refs. [39, 40] are used to simulate the decay dynamics. We refer to these as “non-resonant” $D^{**}$ decays, in contrast to “resonant” $D^{**}$ decays, which proceed via well-defined intermediate states with Breit-Wigner resonance shapes characterized by specific mass and width parameters. Non-resonant $D^{**}$ decays account for approximately $0.7\%$ of all semileptonic and semitauonic events in our reconstructed sample. We also use simulated samples of the isospin-conjugate modes to understand the contamination from $B^{+}$ decays.

The simulation is corrected using data-driven weights to account for differences in identification and reconstruction efficiencies. Lepton identification (LID) efficiency and fake rate corrections for electrons are applied as functions of the laboratory-frame momentum, angle relative to the electron beam, and charge of the electron candidate. These corrections are derived from samples of $e^{+}e^{-}\to e^{+}e^{-}(\gamma)$, $e^{+}e^{-}\to e^{+}e^{-}e^{+}e^{-}$, and events with $J/\psi\to e^{+}e^{-}$ decays. Muon LID corrections are obtained using samples of $e^{+}e^{-}\to\mu^{+}\mu^{-}\gamma$, $e^{+}e^{-}\to e^{+}e^{-}\mu^{+}\mu^{-}$, and events with $J/\psi\to\mu^{+}\mu^{-}$ decays. The rates of misidentifying charged hadrons as leptons are corrected using samples of $K_{S}^{0}\to\pi^{+}\pi^{-}$, $D^{*+}\to D^{0}\pi^{+}$, and $e^{+}e^{-}\to\tau^{+}\tau^{-}$. The efficiency for identifying slow pions from $D^{*+}\to D^{0}\pi^{+}$ decays is corrected using studies of $\overline{B}{}^{0}\to D^{*+}\pi^{-}$. All data and simulated events are reconstructed and analyzed with the open-source basf2 framework [42].

To select events likely to contain $\Upsilon(4S)\to B^{0}\overline{B}{}^{0}$ decays, we require at least three tracks and three ECL clusters in the event. Tracks must have transverse momenta greater than $0.1$ GeV and originate within $|d_{0}|<0.5$ cm and $|z_{0}|<2.0$ cm. Here, $|d_{0}|$ and $|z_{0}|$ denote the distance of closest approach between the nominal interaction point (IP) and the track in the plane perpendicular to and along the beam axis, respectively. At this stage, all tracks are assigned a pion mass hypothesis. We reconstruct ECL clusters with energy deposits above $0.1$ GeV that are not associated with any track. Finally, the sum of the selected track and cluster energies must exceed $4$ GeV.

We reconstruct $B_{\mathrm{tag}}$ candidates using the Full Event Interpretation (FEI) algorithm [26]. The algorithm constructs $B_{\mathrm{tag}}$ candidates from tracks and clusters, using multivariate classifiers and by using a hierarchical approach. The algorithm is trained to identify semileptonic decays and we use it to reconstruct $B^{0}\to D\ell\nu$ and $B^{0}\to D^{*}\ell\nu$ decay candidates, where the $D$ and $D^{*}$ undergo subsequent hadronic decays. A complete list of decay modes and selection criteria is given in Ref. [43]. Each $B_{\mathrm{tag}}$ candidate is assigned a confidence score by the algorithm, ranging from zero to one. Candidates with a confidence score above $0.1$ are selected. To suppress signal-side semitauonic decays in the $B_{\mathrm{tag}}$ candidates, the lepton momentum in the center-of-mass (c.m.) frame must exceed $1$ GeV. The cosine of the angle between the $B$ meson’s momentum and its visible decay products in the c.m. frame is defined as

where $E_{\mathrm{beam}}$ is the beam energy, $m_{B}$ is the $B$ meson mass, and $|\vec{p}{}_{B}|$ is its momentum, computed from $m_{B}$ and $E_{\mathrm{beam}}$. Here $Y=D^{(*)}\ell$ represents the system of visible decay products. For correctly reconstructed semileptonic $B$ decays with a single undetected neutrino, $\cos\theta_{BY}$ lies within $[-1,1]$, but resolution effects and final-state radiation shift it beyond this range. Semitauonic decays with multiple missing neutrinos have on average large negative values of $\cos\theta_{BY}$. All $B_{\mathrm{tag}}$ candidates are required to have $\cos\theta_{BY}$ in the range $[-1.75,1.1]$, reducing the fraction of semitauonic decays among selected candidates to below $0.4\%$.

Background from $q\bar{q}$ production is suppressed using Fox-Wolfram moments [44], which are constructed from a superposition of spherical harmonics using tracks and clusters. The tracks used in the calculation of Fox–Wolfram moments must lie within the CDC acceptance region, have transverse momenta above $0.1$ GeV, and satisfy $|z_{0}|<3.0$ cm and $|d_{0}|<0.5$ cm. The less stringent $|z_{0}|$ requirement enhances discrimination against $q\bar{q}$ backgrounds by including tracks with larger longitudinal displacements, which are more characteristic of signal decays and improve event shape information. We apply a cut on the ratio of the second-order to the zeroth-order Fox–Wolfram moments, with the second-order moment quantifying deviations of the energy flow from isotropy and the zeroth-order moment reflecting the event’s spherical geometry. Higher ratios indicate a more collimated structure typical of $q\bar{q}$ events. Therefore, we require this ratio to be less than 0.4 to reduce the $q\bar{q}$ background.

On average, approximately 2.02 $B_{\mathrm{tag}}$ candidates per event are reconstructed in the events that passed the selection criteria. From these, we select the candidate with the highest classifier value.

We reconstruct $B_{\mathrm{sig}}$ candidates in two final states: $D^{+}\ell^{-}$ and $D^{*\,+}\ell^{-}$ candidates are formed from tracks and clusters not associated with the $B_{\mathrm{tag}}$ candidate. Exactly one charged lepton candidate is required. The signal-side lepton is required to have a charge opposite to the $B_{\mathrm{tag}}$ lepton and is identified using a likelihood-based score that incorporates information from several subdetectors. Electron identification relies on information from the ECL, CDC, TOP, and ARICH, with the most important discriminant being the ratio of reconstructed ECL energy to the estimated track momentum, which is expected to be close to unity for electrons. Electron candidates must have momenta above $0.2$ GeV, with a loose LID score selection. We correct the electron energy for bremsstrahlung losses by adding back ECL clusters that are near the tracks, following the methodology of Ref. [24]. Muons are identified by extrapolating tracks to the KLM, where the likelihood is primarily constructed from the longitudinal penetration depth and transverse scattering of the extrapolated track. Muon candidates must have momenta above $0.4$ GeV, with a stringent LID score requirement. The efficiency for correctly identifying electrons is 98.7% (99.7%) in semitauonic (semileptonic) events, with a misidentification rate such that only 1% of pions or kaons pass this requirement. The efficiency for correctly identifying muons is 79.3% (86.3%) in semitauonic (semileptonic) events, with only 5% of pions or kaons satisfying the selection.

Neutral pion ($\pi^{0}$) candidates are reconstructed from pairs of ECL clusters not associated with any tracks with an invariant mass between $120$ and $145$ MeV. To suppress background, clusters must have energies above $0.08$, $0.03$, or $0.06$ GeV in the forward, barrel, and backward regions, respectively. Each cluster must consist of multiple crystals, lie within the CDC angular acceptance, and have a measured time within $200$ ns of the expected event time. A multivariate classifier is constructed from electromagnetic cluster shape quantities and combined with the photon’s distance to the nearest track to distinguish real photons from clusters originating from hadronic showers. More details can be found in Ref. [23].

Neutral kaon ($K_{S}^{0}$) candidates are reconstructed from pairs of charged particles, each of which is assigned a pion mass hypothesis, whose combined invariant mass lies between $450$ and $550$ MeV and which can be fit to a common vertex. The flight distance must be positive, the significance of the displacement between the point of closest approach and the IP must exceed $0.5$, and the cosine of the angle between the momentum and vertex position vector must be greater than $0.8$. The displacement significance, defined as the separation divided by its uncertainty, distinguishes true $K_{S}^{0}$ candidates from random combinations of tracks.

The decays of $D$ mesons are reconstructed in modes with large branching fractions and high purities. For the $D^{+}\ell^{-}$ final state, we include $D^{+}\to K^{-}\pi^{+}\pi^{+}$, $K^{0}_{S}\pi^{+}\pi^{0}$, $K^{0}_{S}\pi^{+}\pi^{+}\pi^{-}$, $K^{0}_{S}\pi^{+}$, $K^{-}K^{+}\pi^{+}$, and $K^{0}_{S}K^{+}$. For the $D^{*+}\ell^{-}$ final state, where only the $D^{*+}\to D^{0}\pi^{+}$ decay is used, we include $D^{0}\to K^{-}\pi^{+}\pi^{0}$, $K^{-}\pi^{+}\pi^{+}\pi^{-}$, $K^{0}_{S}K^{+}K^{-}$, $K^{+}K^{-}$, $K^{-}\pi^{+}$, $K^{0}_{S}\pi^{+}\pi^{-}$, and $\pi^{-}\pi^{+}$. All $D$ candidates must have a reconstructed invariant mass within $2.5\sigma$ of the nominal $D$ mass [41], where $\sigma$ refers to the resolution of the mass peak.

For the reconstruction of $D^{*+}$ candidates, each $D^{0}$ candidate is combined with a single charged track, assumed to be a pion. Candidates are required to have a mass difference $\Delta m=m(D^{*+})-m(D^{0})$ between 130 and 160 MeV, corresponding to 2.8 times the $\Delta m$ resolution. An additional vertex fit is performed on each $D^{*+}$ candidate to update its momentum. Using MC simulation, we estimate that for true $D^{*+}$ candidates, the slow pion is correctly identified in 71% of cases.

The construction of $B_{\mathrm{sig}}$ candidates is achieved by combining $D^{(*)+}$ and $\ell^{-}$ candidates. The requirement $-15<\cos\theta_{BY}<1.1$ selects semitauonic and semileptonic final states, with approximately 5% of semitauonic signal events falling outside this range on the negative end.

Finally, $B_{\mathrm{tag}}$ and $B_{\mathrm{sig}}$ candidates are combined to form $\Upsilon(4S)$ candidates, requiring no additional charged tracks in the event. We demand that $B_{\mathrm{tag}}$ and $B_{\mathrm{sig}}$ are of opposite flavor and reject events in which one $B$ mixed, as such events possess lower purity. On average, fewer than $1.08$ $\Upsilon(4S)$ candidates remain per event after applying the selection criteria, and we select a single candidate per event based on the highest signal-side $D^{(*)+}$ vertex fit $p$-value. The unassigned energy in the calorimeter of the selected $\Upsilon(4S)$ candidate, $E_{\mathrm{extra}}$, is calculated by summing clusters not associated with any particles used in the reconstruction. Two multivariate algorithms, described in Ref. [45] and based on the FastBDT classifier [46], remove contributions from beam background and hadronic split-offs, by utilizing features based on timing, energy spread, scintillation pulse shape, and cluster localization in the ECL. Clusters added as bremsstrahlung corrections to tracks are excluded. From correctly reconstructed events we expect values of $E_{\mathrm{extra}}$ near zero, while background events typically exhibit higher values. We only retain events with $E_{\mathrm{extra}}<1.2$ GeV.

Appendix A provides more details on the differences of the signal and tag-side selection and reconstruction.

Signal extraction is performed using a multiclass classification algorithm that differentiates semitauonic signal, semileptonic signal, and background events. The model employs gradient-boosted decision trees (BDTs), where each tree corrects the errors of the previous ones to improve classification. Further details can be found in Ref. [47].

The BDT is trained on five input variables. The most discriminating variable is $\cos\theta_{BY}$ of the $B_{\mathrm{sig}}$ candidate, followed by the unassigned energy in the calorimeter, $E_{\mathrm{extra}}$. The third most important variable, $\cos^{2}\Phi_{B}$, is defined as

It combines $\cos\theta_{BY}$ from both the $B_{\mathrm{sig}}$ and $B_{\mathrm{tag}}$ candidates with the angle $\gamma$ between their $Y$ momenta. For correctly reconstructed semileptonic $B_{\mathrm{sig}}$ and $B_{\mathrm{tag}}$ candidates, each with a single missing neutrino, $\cos^{2}\Phi_{B}$ is expected to take values between zero and one. Events with multiple missing neutrinos, such as semitauonic decays or misreconstructed events, tend to have larger values. The fourth and fifth most important input variables are the center-of-mass momenta of the $D$ ($p_{D}^{*}$) and lepton ($p_{\ell}^{*}$) candidates, respectively. These variables help distinguish semitauonic, semileptonic, and background events based on the different phase space available in each case. Figure 1 shows the five input variables for $D$ and $D^{*}$ candidates with electrons and muons combined. A good separation of all three event types can be obtained between semileptonic and semitauonic signal decays in $\cos\theta_{BY}$, $\cos^{2}\Phi_{B}$, $p_{\ell}^{*}$, and $p_{D}^{*}$, whereas $E_{\mathrm{extra}}$ is very powerful at separating semileptonic and semitauonic signal from other semileptonic processes or backgrounds. The three resulting classification scores are denoted as $z_{\tau}$, $z_{\ell}$, and $z_{\mathrm{bkg}}$ for semitauonic, semileptonic, and background events, respectively.

We extract the signal using a binned two-dimensional log-likelihood fit to the variables $z_{\tau}$ and $z_{\mathrm{diff}}=z_{\ell}-z_{\mathrm{bkg}}$. We consider four categories of events: $D^{+}e^{-}$, $D^{*+}e^{-}$, $D^{+}\mu^{-}$, and $D^{*+}\mu^{-}$ candidates. The likelihood is implemented using the pyhf package [48, 49]. The total likelihood function has the form

with the individual category likelihoods $\mathcal{L}_{c}$ and nuisance-parameter (NP) constraints $\mathcal{G}_{l}$. The product in Eq. 4 runs over all categories $c$ and independent uncertainty sources $l$, respectively. The role of the NP constraints is detailed in Section VII. Each category likelihood $\mathcal{L}_{c}$ is defined as the product of individual Poisson distributions $\mathcal{P}$,

with $n_{i}$ denoting the number of observed data events and $\nu_{i}$ the total number of expected events in a given bin $i$. The number of expected events in a given bin, $\nu_{i}$, is estimated using simulated events. It is given by

where $\eta_{k}$ is the total number of events from a given process $k$ with a fraction $f_{ik}$ of such events being reconstructed in the bin $i$. The values of $\cal{R}(D^{+})$ and $\cal{R}(D^{*+})$ and the sum of semileptonic signal decays $\eta_{D^{(*)}\ell}$ from electrons and muons are used to determine the number of semitauonic decays $\eta_{D^{(*)}\tau}$ via

where $\epsilon_{D^{(*)}\tau}$ and $\epsilon_{D^{(*)}\ell}$ denote the efficiencies of semitauonic and semileptonic signal decays from electrons and muons.

The factor of $\frac{1}{2}$ in Eq. 6 accounts for the definition of the semileptonic signal, which includes both electron and muon final states. We implement a non-uniform binning scheme in both dimensions, increasing granularity in regions where the classifier is most sensitive to semitauonic events and $\overline{B}{}\to D^{**}\ell\bar{\nu}_{\ell}$ as well as other semileptonic gap backgrounds. In contrast, the region dominated by semileptonic $\overline{B}{}\to D^{(*)}\ell\bar{\nu}_{\ell}$ events is binned more coarsely, as these events are readily identifiable. Figure 2 shows the binning structure for semitauonic and semileptonic events, for the background category from $\overline{B}{}\to D^{**}\ell\bar{\nu}_{\ell}$ and gap modes, and for other backgrounds.

The likelihood (Eq.4) is numerically maximized to fit the values of the different components, $\eta_{k}$, using the observed events. This maximization is performed with the iMinuit package[50]. The ten parameters of interest we determine are:

• •

  $\mathcal{R}(D^{+})$ and $\mathcal{R}(D^{*+})$ (2 parameters, shared between $e$ and $\mu$ channels);

• •

  Normalizations $\eta_{D^{(*)}\ell}$ of $\overline{B}{}^{0}\to D^{+}\ell^{-}\bar{\nu}_{\ell}$ and $\overline{B}{}^{0}\to D^{*+}\ell^{-}\bar{\nu}_{\ell}$ events (2 parameters, shared between $e$ and $\mu$ channels);

• •

  Background from $B\to D^{**}\ell\bar{\nu}_{\ell}$ and semileptonic gap events (2 parameters; shared between $e$ and $\mu$ channels);

• •

  Number of other $B\overline{B}{}$ or continuum background events (4 parameters; one for each fit category).

We also test LFU between electrons and muons for the semileptonic signal modes. For this we modify the fit setup by introducing four normalization parameters for each of the four decays of interest and modify Eq. 7 accordingly.

To validate the fit procedure, we generated ensembles of pseudoexperiments for different input $\mathcal{R}(D^{+})$ and $\mathcal{R}(D^{*+})$ values. Fits to these ensembles show no biases in central values and no under- or overcoverage of determined confidence intervals. Assuming SM values for $\mathcal{R}(D^{+})$ and $\mathcal{R}(D^{*+})$, we expect significances of 3.9 and 7.0 standard deviations, respectively.

Several systematic uncertainties affect the measured ratios of $\mathcal{R}(D^{+})$ and $\mathcal{R}(D^{*+})$, and Table 1 provides a summary. The effect of systematic uncertainties is directly incorporated into the likelihood function. We distinguish between additive and multiplicative uncertainties: Additive uncertainties affect the signal and background template shapes, whereas multiplicative uncertainties affect efficiencies or branching fractions of semitauonic and semileptonic signal events. A vector of NPs, $\bm{\theta}$, is introduced and each NP is constrained in the likelihood Eq. 4 using a standard normal distribution $\mathcal{G}_{l}=\mathcal{G}(\theta_{l})$ with $l$ denoting an independent uncertainty source or if applicable Poisson constraints.

A brief summary of each significant uncertainty source follows, ordered by their importance:

The dominant uncertainty arises due to the finite size of the simulated data samples and is estimated using the “Barlow-Beeston Lite” method [49], introducing Poisson constraints to correctly treat sparsely populated bins.

The next largest uncertainty stems from the limited knowledge on the composition and modeling of the semileptonic gap processes. To account for the former, we assign a 100% uncertainty to their branching fractions. The latter is addressed by varying the heavy-quark-symmetry-based form factors independently for each assumed gap process, using the uncertainties and correlations provided in Refs. [39, 40] for the broad states. Due to the larger contamination, their impact is more pronounced in the $D$ channel than in the $D^{*}$ channel, translating into a larger systematic uncertainty for $\mathcal{R}(D^{+})$.

In contrast, the branching fractions and modeling of resonant $\overline{B}{}\to D^{**}\ell\bar{\nu}_{\ell}$ decays are better constrained experimentally [51] and we vary the form factors of the broad and narrow states using the uncertainties and correlations provided in Refs. [39, 40].

Lepton identification impacts the precision of $\mathcal{R}(D^{+})$ and $\mathcal{R}(D^{*+})$ due to the limited size and systematic uncertainties of the calibration samples used to correct discrepancies between data and Monte Carlo simulations. These effects influence both correctly identified leptons and misidentified leptons (“fakes”). The contribution of misidentified background events is larger in the signal-enriched regions of the $D(e/\mu)$ channels compared to the $D^{*}(e/\mu)$ channels, resulting in a higher uncertainty in $\mathcal{R}(D^{+})$.

We assign a track reconstruction efficiency uncertainty of 0.3% per track for kaon, pion, and lepton tracks, reflecting the imperfect knowledge of the track-reconstruction efficiency. This uncertainty is estimated using a control sample of $e^{+}e^{-}\to\tau^{+}\tau^{-}$ events and is assumed to be fully correlated across all tracks. The slow-pion efficiency for momenta less than 0.2 GeV is corrected relative to the tracking efficiency at momenta larger than 0.2 GeV in the laboratory frame. We study $\overline{B}{}^{0}\to D^{*+}\pi^{-}$ decays to determine corrections for three momentum bins spanning $[0.05,0.12,0.16,0.20]$ GeV. The associated statistical and systematic uncertainties of the correction weights affect the $\mathcal{R}(D^{+})$ and $\mathcal{R}(D^{*+})$ and are evaluated using variations of the correction weights, taking into account their correlations.

A systematic uncertainty also arises from the modeling of the BDT input variables. The most significant mis-modeling is observed in the negative region of the $\cos\theta_{BY}$ distribution of the $B_{\mathrm{sig}}$ candidate. To estimate the impact of this mis-modeling, MC samples are reweighted to data, with event weights derived from a linear spline fit to the data-to-MC ratio as a function of $\cos\theta_{BY}$. This reweighting is done separately for each of the four categories, yielding distinct spline fits. The event weights are then used to construct a fit, with each template normalized to the nominal fit event count to isolate shape effects. Two approaches were investigated: first, reweighting data while fitting with nominal templates, and second, reweighting the fit templates while fitting to nominal data. We use the larger shift to assess the uncertainty.

The uncertainty on the form factors of semitauonic and semileptonic signal are evaluated by using the eigenvariations provided in Ref. [6]. In addition, we assign the difference in central value between the parametrization of Ref. [6] and Refs. [52, 53] for semileptonic signal and Ref. [54] for semitauonic signal events, using the information from Refs. [55, 56].

Finally, we assign a 10% uncertainty to the fraction of continuum events in the background, based on the observed difference in efficiency between off-resonance data and the expected continuum contribution.