Existing measurements of jV cb j use either exclusive final states with B → D Ã lν l and B → Dlν l providing the most precise values or inclusive final states. In inclusive determinations of jV cb j, the total decay rate can be expressed as an expansion of a small number of nonperturbative matrix elements with the heavy-quark expansion (HQE). Using HQE, the total semileptonic rate can be expanded in powers of Λ QCD =m b , the ratio of the QCD scale parameter and the bottom-quark mass and perturbative corrections proportional to the strong coupling constant α s can also be systematically incorporated [1][2][3][4][5][6][7][8].
The current world averages [9] for jV cb j determined from inclusive and exclusive approaches are jV incl cb j ¼ ð42.19 AE 0.78Þ × 10 -3 and ð1Þ
jV excl cb j ¼ ð39.25 AE 0.56Þ × 10 -3 ; ð2Þ
respectively. The uncertainties are the sum of experimental and theoretical uncertainties; the world averages differ by about 3 standard deviations. The 2% relative uncertainty in the world average for the inclusive approach is largely due to the theory uncertainty associated with the truncation of HQE and perturbative expansion [10,11]. To further reduce this uncertainty, higher order nonperturbative matrix elements must be determined from measured spectral moments. This is complicated by the proliferation of HQE parameters at higher orders in the expansion. At Oð1=m 4 b Þ in the HQE 13 nonperturbative matrix elements contribute to the total rate and the spectral energy and mass moments.
Reference [12] outlines a novel and alternative approach to determine jV cb j from inclusive decays avoiding this proliferation of terms. Exploiting reparametrization invariance, the authors reduce the number of parameters necessary to calculate the total rate at Oð1=m 4 b Þ to only eight. Unfortunately, spectral moments of lepton-energy and hadron-mass spectra violate reparametrization invariance. However, reparametrization invariance is retained in the spectral moments of the lepton mass squared 2 where p i is the four-momentum of the particle i.
We present measurements of the spectral moments of the lepton mass squared hq 2n i with n ¼ 1 -4 for q 2 > 1.5 GeV 2 =c 4 up to 8.5 GeV 2 =c 4 . The simultaneous analysis of these moments can determine the nonperturbative matrix elements as their contributions vary with the q 2 threshold [12]; moments with a lower q 2 threshold retain more information about the inclusive B → X c lν l process. Charge conjugation is implied throughout this paper, and BðB → X c lν l Þ is defined as the average of the branching fraction with B 0 and B þ and l ¼ e, μ.
We present raw and central moments, with the latter having the benefit of smaller correlations between q 2 thresholds and the orders of moments. The first measurement of the first q 2 moment was reported in Ref. [13] with an implicit lower requirement on the lepton energy of 1 GeV. This requirement renders the measured moment unsuitable for the analysis outlined in Ref. [12].
A measurement of the q 2 moments, similar to the one presented in this paper, using the full Belle data set was recently reported by the Belle Collaboration [14] for q 2 > 3.0 GeV 2 =c 4 , covering 58% of the accessible B → X c lν l phase space. We report measurements of the raw and central q 2 moments with comparable precision and include for the first time the experimentally challenging low q 2 region q 2 > 1.5 GeV 2 =c 4 , covering up to 77% of the accessible B → X c lν l phase space.
The remainder of this paper is organized as follows: Section II describes the data set used in this analysis, the Belle II 2 detector, and the simulation of e þ e -collision events. Section III introduces the tag-side and the inclusive reconstruction of semileptonic B decays. Section IV describes the background subtraction, calibration, and calculation of the lepton mass squared moments. Section V discusses the systematic uncertainties affecting the measurement. Section VI presents the main findings, and Sec. VII contains our conclusions.
We analyze data collected in 2019 and 2020 by the Belle II 2 detector [15] at the SuperKEKB e þ e -accelerator complex [16]. At SuperKEKB, 7 GeVelectrons collide with 4 GeV positrons giving a c.m. energy of ffiffi ffi s p ¼ 10.58 GeV, corresponding to the mass of the ϒð4SÞ resonance. This results in a boost of βγ ¼ 0.28 of the c.m. frame relative to the laboratory frame. The integrated luminosity of 62.8 fb -1 [17] of the data corresponds to ð68.2 AE 0.9Þ × 10 6 B pairs. We use 9.2 fb -1 of data recorded 60 MeV below the ϒð4SÞ resonance to constrain contributions from e þ e -→ q q continuum processes.
The Belle II 2 detector is a substantial upgrade of the Belle detector [18] with improved reconstruction of charged and neutral particles and particle identification performance. The detector consists of several subdetectors arranged in a cylindrical structure around the e þ e -interaction point (IP). The IP is enclosed by a beryllium beam pipe with an inner radius of 1 cm. The part of the detector closest to the IP is the pixel detector (PXD), consisting of two layers of depleted p-channel field-effect-transistor pixel-sensor modules [19]. The first layer comprises sixteen modules arranged in eight ladders. The second layer was only partially installed for data taking and consists of four modules. The PXD is surrounded by four layers of doublesided silicon strip modules: the silicon vertex detector (SVD). The first SVD layer is arranged parallel to the beam axis, while the forward sections of the second to fourth layers are tilted with respect to the beam axis in order to reduce the overall material budget and the number of sensors. Both silicon tracking detectors are enclosed by the central drift chamber (CDC), which is filled with a He (50%) and C 2 H 6 (50%) gas mixture. The CDC contains 56,576 sense and field wires oriented along the beam direction or tilted and arranged into 56 radial layers.
By combining the information from axial and stereo wires, the full three-dimensional trajectory of a charged particle is reconstructed, and its specific ionization dE=dx is measured. Outside the CDC, a time-of-propagation detector (TOP) and an aerogel ring-imaging Cherenkov detector (ARICH) cover the barrel and forward endcap regions of the detector, respectively. The TOP reconstructs spatial and temporal coordinates of the ring of Cherenkov light cones emitted from charged particles passing through quartz radiator bars. The information from both the TOP and ARICH and the CDC are combined together to identify charged particles. The electromagnetic calorimeter (ECL) consists of a 3 m long barrel section with an inner radius of 1.25 m and annular endcaps. In total 8736 CsI(Tl) crystals arranged in a pointing geometry allow for precise energy and timing measurements of neutral and charged particles. The ECL is located outside the TOP and inside the remaining volume of a superconducting solenoid with a field strength of 1.5 T. The K 0 L and muon detector (KLM) is located outside of the coil. It consists of an alternating structure of 4.7 cm thick iron plates and active detector elements. The iron plates are used as the magnetic flux return yoke for the solenoid and absorber material to range out charged hadrons. The detector elements are glasselectrode resistive plate chambers and plastic scintillators in the barrel and endcap regions, respectively.
We define the z axis of the laboratory frame as the central axis of the solenoid with the positive direction in the direction of the electron beam. The polar angle θ and the longitudinal and transverse directions are defined with respect to the z axis. Variables with asterisk superscripts are measured in the c.m. frame; variables without asterisks are measured in the laboratory frame.
Charged particle tracks are reconstructed by combining information from the PXD, SVD, and CDC [20]. The reconstruction of energy depositions from neutral and charged particles in the ECL (ECL clusters) uses shower shapes and timing information [21]. Tracks are identified as electron or muon candidates by combining information from several subdetectors into a single lepton identification likelihood L l (PID). Muons are identified reliably by extrapolating tracks to the KLM. The main features used for the construction of the likelihood are the longitudinal penetration depth and the transverse scattering of the extrapolated track in the KLM. For electrons, the likelihood is constructed from information from the ECL, CDC, TOP, and ARICH. The most important discriminant is the ratio of the reconstructed energy in the ECL to the estimated track momentum, which should be close to unity for electrons. The identification of charged pions, kaons, and protons is based on likelihood information from the CDC, TOP, and ARICH. Their likelihoods are denoted as L π , L K , and L p . Hadrons with momenta less than 700 MeV=c are primarily identified using dE=dx measurements from the CDC. Hadrons with momenta larger than 700 MeV=c are primarily identified using the TOP and ARICH measurements. Photon candidates are identified using the ECL shower shape of clusters not matched to a track. We require each photon candidate to have a transverse energy greater than 30 MeV when reconstructed in the barrel or 20 MeV when reconstructed in either endcap. A loose selection on a multivariate shower-shape classifier that uses multiple Zernike moments [22] is imposed. A more detailed overview of the Belle II 2 PID algorithms and the photon reconstruction algorithms can be found in Ref. [21].
Monte Carlo (MC) samples are used to determine reconstruction efficiencies and acceptance effects as well as to estimate background contamination. MC samples of B decays are simulated using the EvtGen generator [23]. The simulation of e þ e -→ q q continuum processes is carried out with KKMC [24] and PYTHIA8 [25]. Electromagnetic final-state radiation (FSR) is simulated using PHOTOS [26] for all charged final-state particles. Interactions of particles with the detector are simulated using GEANT4 [27].
The simulation is corrected using data-driven weights to account for differences in identification and reconstruction efficiencies. The PID for electrons is corrected as a function of the laboratory-frame momentum and polar angle and charge of the electron candidate using samples of e þ e -→ e þ e -ðγÞ and e þ e -→ e þ e -e þ e -events and events with J=ψ → e þ e -decays. The PID for muons is corrected using samples of e þ e -→ μ þ μ -γ and e þ e -→ e þ e -μ þ μ -, and events with J=ψ → μ þ μ -decays. The average multiplicative corrections are 0.95 and 0.89 for electron and muon candidates, respectively. The rates of misidentifying charged hadrons as charged leptons are corrected using samples of
, and e þ e -→ τ þ τ -, with average multiplicative misidentification-rate corrections of 1.50 and 0.98 for electron and muon candidates, respectively.
All recorded e þ e -collision data and simulated events are reconstructed and analyzed with the open-source basf2 framework [28].
The analysis relies on accurate modeling of B → X c lν l decays. Inclusive semileptonic B → X c lν l decays are dominantly B → Dlν l and B → D Ã lν l decays. The B → Dlν l decays are modeled using the BGL parametrization [29] with form-factor parameter values and uncertainties from the fit in Ref. [30]. For B → D Ã lν l decays, the BGL implementation proposed in Refs. [31,32] with form-factor parameter values and uncertainties from a fit to the measurement of Ref. [33] is used. Both branching fractions are normalized to the average branching fraction of Ref. [9] assuming isospin symmetry. Semileptonic B → D ÃÃ lν l decays with D ÃÃ ¼ D Ã 0 ; D 0 1 ; D 1 ; D Ã 2 are modeled using heavy-quark-symmetrybased form factors proposed in Ref. [34] and with D ÃÃ masses and widths from Ref. [35].
For the B → D Ãà lν l branching fractions, we adopt the values of Ref. [9] to account for missing isospin-conjugated and other established decay modes observed in studies of B decays into fully hadronic final states. This follows the prescription outlined in Ref. [34]. All existing exclusive B → D Ãà lν l measurements only use D ÃÃ0 → D ðÃÞþ π - decay modes. To correct for the missing isospin modes we multiply the branching fractions with a multiplicative factor of 3=2.
In the average in Ref. [9], all measurements of B → D Ã 2 lν l are relative to DÃ 2 → D Ã-π þ . To account for DÃ 2 → D -π þ contributions, we apply a multiplicative factor of 1.54 AE 0.15 calculated from the branching fractions of Ref. [35].
The world average for B → D 0 1 lν l in Ref. [9] combines measurements that only marginally agree with each other (the probability of the combination is below 0.01%). We exclude the measurement of Ref. [36] that is in conflict with the measured branching fractions of Refs. [37,38]. That measurement also conflicts with the expectation that 39,40]. By excluding Ref. [36] we obtain
The world average for BðB → D 1 lν l Þ does not include contributions from D 1 → Dππ. To account for these, we use a multiplicative factor 0.43 AE 0.11 calculated from the branching fractions of [41]. The contribution of D 1 → Dππ decays is subtracted from the B → Dππlν l branching fraction measured in Ref. [42]. The three-hadron final states must be corrected for missing isospin-conjugated modes. Following Ref. [42], we use an average isospin correction multiplicative factor of
whose uncertainty covers the isospin hypotheses for different resonant final states [f 0 ð500Þ → ππ and ρ → ππ result in f ππ ¼ 2=3 and 1=3, respectively] and nonresonant threebody decays (f ππ ¼ 3=7). Furthermore, it is assumed that the resulting branching fractions saturate the branching fractions of orbitally excited states:
For the B → D ðÃÞ ππlν l contributions not covered by decays into D 1 → Dππ, we use values measured in Ref. [42]. We neglect the small contribution from B → D ðÃÞ s Klν l decays. There is still a difference between the sum of all exclusive modes and the inclusive B → X c lν l branching fraction of Ref. [35]. In the following, this missing component contributing to the total branching fraction is referred to as the "gap." We fill this gap with equal parts of B → Dηlν l and B → D à ηlν l decays and assign an uncertainty of 100% to its branching fraction. These decays are simulated with final-state momenta uniformly distributed in the available phase space or an alternative model involving a broad resonance for the hadronic X c final state.
Figure 1 shows the resulting q 2 spectrum evaluated without reconstruction effects for the different X c final states, and Table I summarizes the semileptonic branching fractions. At high q 2 , contributions from B → D à lν l dominate, whereas at low q 2 , B → D Ãà lν l and nonresonant X c (B → D ðÃÞ ππlν l and gap processes) have sizable contributions.
A. Tag-side reconstruction
We reconstruct ϒð4SÞ → B B events with the full event interpretation (FEI) algorithm [43]. The algorithm reconstructs one of the B mesons of the B B pair in fully hadronic decays. In the following, the tag-side B candidate reconstructed by the FEI is denoted as B tag . The FEI uses a hierarchical bottom-up approach starting with the selection of charged and neutral final-state particles (e -, μ -, π -, K -, p, γ) from tracks, and ECL clusters, combining them into intermediate particles
), and finally forming B tag candidates. At each stage, the FEI uses an optimized implementation of gradient-boosted decision trees [44] to estimate the signal probability P FEI of each candidate in a distinct decay chain to be correctly reconstructed. For each candidate, the decision trees combine the signal probability of previous stages with additional kinematic and vertex-fit information. More than 100 decay channels are reconstructed, resulting in Oð10;000Þ decay chains.
We select events that have at least three charged particles and three ECL clusters to suppress B tag candidates from continuum processes. The total visible energy of the event in the c.m. frame must be greater than 4 GeV, and the total energy in the ECL is required to be between 2 and 7 GeV. To reduce continuum background, events must have R 2 < 0.4, with R 2 the ratio of the second and zeroth Fox-Wolfram moments [45]. We suppress continuum events by requiring cosðθ T Þ < 0.7, where θ T is the angle between the thrust axis of the decay products of the B tag and the thrust axis of the rest of the event [46]. Note that B tag candidates are selected by requiring P FEI > 0.01. The reconstruction efficiencies with this requirement are approximatively 0.26% and 0.35% for neutral and charged B tag candidates, respectively. More details on the FEI performance with Belle II 2 data can be found in Ref. [47].
We require B tag candidates to have beam-constrained mass values satisfying
The q 2 spectrum for different X c final states without reconstruction effects ("gen"). Details about the simulation are given in the text. where p à B tag is the three-momentum of the B tag candidate. The energy difference
must be within ½-0.15; 0.1 GeV, where E Ã B tag is the energy of the B tag . All tracks and ECL clusters not used in the reconstruction of the B tag candidate are used to define and reconstruct the signal side. At this stage, we allow for multiple B tag candidates in each event.
Semileptonic B decays are identified by selecting electron and muon candidates with laboratory frame momenta greater than 0.5 GeV=c. These tracks are required to originate from the IP by requiring d r < 1 cm and jd z j < 2 cm. Here, d r and d z are the distances of closest approach to the IP transverse to and along the z axis, respectively. Each lepton candidate is required to have a polar angle within the CDC acceptance ½17°; 150° and at least one hit in the CDC.
The momentum and polar angle selection affects the selection efficiency as a function of q 2 , which is illustrated in Fig. 2. At low q 2 thresholds, the efficiency depends on the final states. A lower selection efficiency is observed for the D ÃÃ and nonresonant contributions, introducing a dependence of the moments on modeling of B → X c lν l . To minimize extrapolation of the moments to unmeasured phase-space regions, we require q 2 > 1.5 GeV 2 =c 4 .
Lepton candidates are selected using
and we require P l > 0.9 for both electrons and muons. To account for the energy of electrons lost to bremsstrahlung photons, the four-momenta of such photons are added to the fourmomenta of electrons. Bremsstrahlung photons are identified using the electron track, extrapolating its PXD and SVD hits and the estimated track intersections with the beam pipe and inner wall of the CDC to the ECL to search for clusters. ECL clusters with energies between 2% and 100% of the electron energy and without any other track association are identified as potential bremsstrahlung photons. All clusters that lie within 3 times the expected resolutions in polar and azimuthal angles are used to correct the electron candidate. These clusters are then removed from consideration for the remainder of the analysis. For charged B tag candidates, we require the signal-side lepton to have a charge opposite to that of the B tag .
Particles with transverse momenta less than 275 MeV=c have radii of curvature in the magnetic field sufficiently small that they loop within the CDC volume and may be reconstructed as multiple tracks. To identify such tracks, we compare the proximity and the magnitude of the momenta of all low-momentum tracks. When there are potential duplicates, we select the track with the smallest value of ð5 × d r Þ 2 þ jd z j 2 . The size of the scaling factor on d r is optimized to minimize track duplicates.
After reconstructing the B tag and signal-side lepton candidate, the X c system is identified as the remaining charged particles and photons. The four-momentum for a charged particle is calculated from the reconstructed track momentum and the assigned mass hypothesis based on the largest identification probability. As we do not explicitly reconstruct charmed states, we denote the reconstructed system in the following as X and its four-momentum p X and mass M X . A signal-side candidate is rejected if the X system does not contain at least one charged particle and the absolute event charge is > 1.
The missing four-momentum in the event is reconstructed as
where p e þ e -is the four-momentum of the colliding electron-positron pair. We require E miss > 0.5 GeV and jp miss j > 0.5 GeV=c to improve the resolution on the mass of the hadronic system. The average multiplicity of B tag l candidates is 1.5 per event. In each event, we retain only the one with the highest lepton momentum. When multiple B tag l candidates share the same lepton, one is chosen randomly.
The lepton mass squared is reconstructed as
To improve the resolution of q 2 reco , we exploit the known kinematics of the e þ e - collision and fit for the four-momenta of B tag , X, l, FIG. 2. Selection efficiencies as functions of q 2 threshold q 2 th . The points for different X c final states and the same lower q 2 threshold are shifted horizontally, and the gray and white bands visually group the same q 2 threshold. and ν l . We construct a χ 2 function for each candidate of the form
where pi is the fitted four-momentum, and C i is the covariance matrix of the four-momentum of a given final-state particle. Note that C l is given by the track fit result, while C B tag and C X are estimated using the corresponding four-momentum residuals. Overall, we fit 14 parameters: the four-momenta components of the B tag and X candidates and the momenta components of the signal lepton and neutrino. The energies of the lepton and neutrino are calculated from the momenta assuming p 2 l ¼ m 2 l and p 2 ν ¼ 0. The kinematic fit is then performed by imposing the following constraints,
and
using Lagrange multipliers. For each event the χ 2 function is numerically minimized with the constraints, following the algorithm described in Ref. [48] implemented in SciPy [49].
Figure 3 shows the distribution of the residuals of q 2 before and after the kinematic fit with simulated signal events. Here the residual is calculated from the reconstructed and generated values. The kinematic fit results in more symmetric residuals and a reduction in the tails of the residuals. The rms improves from 5.76 GeV 2 =c 4 to 2.65 GeV 2 =c 4 , and the bias reduces from 3.43 GeV 2 =c 4 to 1.20 GeV 2 =c 4 .
To measure the lepton mass squared moments, background contributions from other processes must be subtracted from the q 2 distribution. Binned likelihood fits are applied to the M X distribution to determine the number of signal and background events. With this information and the shapes of backgrounds from simulation, an event-wise signal probability w is constructed as a function of q 2 reco . Both steps are discussed in Sec. IVA. We correct for acceptance and reconstruction effects by applying an eventwise calibration q 2 reco → q 2 calib and two additional calibration factors C calib and C gen , discussed in Sec. IV B. The background-subtracted q 2 moment of order n is calculated as a weighted mean
with sums over all events. For each reconstructed q 2 threshold, the binned likelihood fit to M X is repeated to update the event-wise signal probability weights. We use thresholds in the range ½1.5; 8.5 GeV 2 =c 4 in steps of 0.5 GeV 2 =c 4 .
The likelihood fit to the binned M X distribution is carried out separately in the B þ l -, B 0 l -, and B 0 l þ channels to account for efficiency differences in the FEI algorithm. Electron and muon channels are not separated. Contributions from B → X u lν l decays are treated as background and have, on average, high q 2 reco . We suppress this background by fitting the range with M X > 0.5 GeV=c 2 . To determine the number of background events in each of these channels as well as for each reconstructed q 2 threshold, we distinguish the following three event categories:
(1) B → X c lν l signal (with yield η sig ), (2) e þ e -→ q q continuum processes (η q q), and
(3) B B background dominated by secondary leptons and hadronic B decays misidentified as signal lepton candidates (η B B). The likelihood is the product of Poisson likelihoods for each bin i with n i observed events and ν i expected events, with
where f ki is the fraction of events of category k reconstructed in bin i as determined with simulated events. The yield η q q is constrained to its expectation as determined from off-resonance data. To reduce the dependence on the FIG. 3. Comparison of reconstructed, fitted, and generated q 2 for B → X c lν l . The residuals are the difference of estimated ("reco") and generated ("gen") values.
modeling of signal and backgrounds, the fit is carried out in five M X bins. For each channel and reconstructed q 2 threshold, an adaptive binning is chosen. The likelihood is numerically maximized using the MINUIT algorithm [50] in scikit-hep/iminuit [51]. The sample composition projections for q 2 reco > 1.5 GeV 2 =c 4 are shown in Appendix A. The M X and q 2 reco distributions with the fitted MC yields are shown in Fig. 4 for q 2 reco > 1.5 GeV 2 =c 4 with finer granularity than used in the fit. The agreement is fair, and the p value from a χ 2 test for the q 2 reco distribution in the range of 1.5 -15 GeV 2 =c 4 is 30%.
The event-wise signal probability w is obtained by constructing a binned probability as a function of q 2 reco via
with n i the observed events in bin i of q 2 reco . Furthermore, fi are the fractions of events for a given background category estimated from the simulation, and η denote the sum of the number of background events from the M X fits.
We calculate a continuous signal probability wðq 2 reco Þ by interpolating the binned distribution with smoothed cubic splines [52]. Negative probabilities are set to zero. The cubic-spline fit and statistical uncertainties of the signal probability are shown in Fig. 5. The statistical uncertainty on hq 2n i is evaluated by a bootstrapping procedure [53], and a selection of spline fits from replicas is shown in Fig. 5. The statistical uncertainty of wðq 2 reco Þ increases towards large q 2 reco .
The q 2 reco distribution is calibrated by exploiting the linear relationship between reconstructed and generated moments. Figure 6 shows the linear relationship for simulated events for the first moment and as functions of q 2 threshold between the reconstructed and true q 2 distribution. We calibrate each event with FIG. 4. The M X and q 2 reco spectra with B → X c lν l and background components normalized to the results of the M X fits are shown for q 2 reco > 1.5 GeV 2 =c 4 . FIG. 5. Binned signal probability w i for q 2 reco > 1.5 GeV 2 =c 4 together with a smoothed cubic-spline fit (dark red). In addition, variations of the signal spline fit (light red) determined with bootstrap replicas are shown.
The first moments are shown as a function of the minimum q 2 requirement on the reconstructed and true underlying q 2 distributions.
with c n and m n the intercept and slope of the linear relationship for a given moment of order n. More details on the linear calibration for the higher moments can be found in Appendix B. Due to the linearity of the calibration, a small bias remains, which we corrected with an additional multiplicative calibration factor in Eq. ( 13) calculated from simulated events by comparing the calibrated hq 2n calib i and true generated hq 2n gen;sel i moments,
The B tag reconstruction and the Belle II 2 detector acceptance and performance result in an additional bias. To account for these effects we apply a second multiplicative calibration factor C gen by comparing the generated moments with all selection criteria applied (hq 2n gen;sel i) to their value without any selection applied (hq 2n gen i),
The hq 2n gen i are determined from a MC sample without PHOTOS simulation which also corrects for FSR.
Both C calib and C gen are determined for each q 2 threshold and from independent samples from those used to determine the linear calibration function. The C calib factors range between 0.98 and 1.02 depending on the reconstructed and generated q 2 threshold. The C gen factors vary between 0.90 and 1.00 with lower q 2 selection threshold values tending to have more sizable corrections. More details on the eventwise calibration can be found in Appendix C.
We use simulated samples to test the robustness of the measurement method and the background subtraction. Closure tests are carried out with ensembles built from independent simulated samples. We observe small deviations of 0.01% to 0.66% caused by imperfections in the interpolation of wðq 2 reco Þ in the extracted q 2 moments. This deviation is treated as a systematic uncertainty; see Sec. V.
We also test the impact of systematically altered generated q 2 shapes for B → X c lν l . The altered shapes are obtained by completely removing the nonresonant B → X c lν l contributions or by applying scaling factors of 2 or 0.5 to the dominant B → Dlν l or B → D Ã lν l contributions. These variations are significantly outside of the quoted uncertainties of Table I. The moments of the samples with the altered generated q 2 shapes are measured with the nominal B → X c lν l composition, and the observed biases are well within the assigned uncertainties.
The consistency of the measurement for electron and muon final states is checked by separately determining the moments; we find good agreement.
Several systematic uncertainties affect the q 2 moments. Their sources can be grouped into two categories. The first consists of systematic uncertainties originating from background subtraction. The fit to the M X distribution assumes the composition of B → X c lν l and relies on data-driven corrections. These and other uncertainties affect wðq 2 reco Þ and must be propagated to the moments. The second category of uncertainties is related to assumptions when calibrating the moments. Modeling of B → X c lν l and of the Belle II 2 detector affects the calibration function and the calibration factors. To assess the effect of each uncertainty source, we derive alternative sets of moments based on either a varied signal probability function or modified calibration. The deviation from the nominal result is used to estimate the systematic uncertainty.
We include uncertainties from the signal and background compositions, MC statistics, and the data-driven correction factors directly into the likelihood function of the M X fit. This is achieved by introducing nuisance parameters θ ki for event category k and bin i, which are constrained with multivariate Gaussians in the likelihood. The fraction of events is replaced in Eq. ( 14) by
and σ ki denotes the uncertainty on the fraction for event category k and bin i.
The composition uncertainties of B → X c lν l are determined with the branching fraction uncertainties listed in Table I. We evaluate the uncertainties of the BGL formfactor parameters for B → Dlν l , B → D Ã lν l using a set of orthogonal parameter variations for each decay. We include the uncertainty of the B → X u lν l branching fraction from Ref. [35]. The efficiencies for identifying or misidentifying leptons and hadrons are estimated from ancillary measurements. We assign a track selection efficiency uncertainty of 0.69% per track on the signal side.
We propagate uncertainties on PID and tracking efficiencies, the B → X u lν l branching fraction, and the background yield obtained from the M X fit to w i ðq 2 reco Þ with all uncertainties varied according to a multivariate Gaussian distribution. We repeat the analysis with varied histograms and take the variation of the resulting moments as the systematic uncertainties due to these sources.
We study the impact of the choice of the smoothing factor for the interpolation of the cubic splines used to derive wðq 2 reco Þ and find it to be negligible.
The calibration curves depend on the composition and modeling of B → X c lν l . We evaluate the impact of the branching fraction uncertainties in B → Dlν l , B → D à lν l , and B → D Ãà lν l by independently varying the branching fraction of each simulated component by 1 standard deviation and determining the corresponding variations of the calibration functions and calibration factors. To assess the effect of the poorly known nonresonant and gap modes, calibration procedures from two different approaches are compared. The first model removes contributions from B → D ðÃÞ ππlν l and B → D ðÃÞ ηlν l decays. The second model replaces them with decays to D Ãà states (D à 0 and D 0 1 ). Although there is no experimental evidence for additional decays of charm 1P states into other final states or the existence of an additional broad state in semileptonic transitions, this provides an alternative kinematic description of the three-body decay, B → D Ãà gap lν l . We also evaluate the sensitivity of the calibration functions and factors to the B → Dlν l and B → D à lν l BGL form-factor parameters. For each orthogonal variation of the BGL parameters we repeat the calibration.
Modeling of the photon and charged-particle multiplicities directly affects the resolution on q 2 and contributes a systematic uncertainty caused by differences between data and MC in how final-state particles are assigned to the signal and tag side. We select a signal-enriched region by requiring M X < 3.0 GeV=c 2 and p à l > 1 GeV=c and calculate correction factors for both multiplicities independently.
We observe differences between data and MC in E miss -jp miss j. We parametrize the differences using a smoothed cubic spline and correct MC events to evaluate the impact on the calibration.
We evaluate the uncertainty from the track finding efficiency and of PID efficiency on the calibration curves.
We propagate the statistical uncertainty on the parameters of the calibration function by varying the calibration curve parameters by 1 standard deviation. For the calibration factors, we vary the statistical uncertainty on C calib × C gen within 1 standard deviation and repeat the calculation of the q 2 moments.
The deviation from the closure for the measurement method discussed in Sec. IV C is assigned as an uncertainty. Its size is subdominant for all moments.
Figure 7 shows the relative systematic uncertainty for the raw moments. A more detailed breakdown of the relative systematic uncertainties is given in Appendix D. For each moment, the total systematic uncertainty decreases with FIG. 7. Total (gray) and grouped (colored histograms) relative systematic uncertainties of the raw q 2 moments as functions of the q 2 threshold. increasing q 2 threshold, whereas the statistical uncertainty increases. At low q 2 thresholds and for the first and second moments, the q 2 reco resolution from mismodeling of the number of charged particles in the X system, the B → X c lν l modeling, and the uncertainty from the background subtraction are of similar size.
The branching fraction and BGL parameter uncertainties of the resonant decays B → Dlν l and B → D Ã lν l are smaller than the uncertainty due to the composition of the higher mass states of the X c spectrum.
At high q 2 thresholds, MC simulation statistics also can be sizable sources of uncertainty for the first and second moments. For the third and fourth moments, the dominant uncertainty at high q 2 thresholds is from the mismodeling of the number of charged particles in the X system, followed by MC simulation statistics and B → X c lν l modeling.
The hq 2n i moments for n ¼ 1 -4 are shown in Fig. 8 for q 2 thresholds ranging from 1.5 GeV 2 =c 4 to 8.5 GeV 2 =c 4 in 0.5 GeV 2 =c 4 increments. Numerical values are given in Appendix D in Tables II-V. Moments with similar q 2 thresholds are strongly correlated. The estimated correlation coefficients are given in Appendix E.
Figure 8 also shows the moments calculated from the simulated B → X c lν l sample constructed with the assumptions described in Sec. II D. The simulated moments include uncertainties from the B → X c lν l composition and B → D ðÃÞ lν l BGL-form-factor parameters. We observe a fair agreement between measured and simulated moments. We compare the raw moments for each order with the simulated moments using χ 2 tests. To obtain numerically stable results, each test only includes measurements with correlation below 95%. The resulting p values range from 27% to 94%.
We calculate values for the central q 2 moments by expanding the binomial relation
and applying the following nonlinear transformation:
FIG. 8. The q 2 moments (blue) as functions of q 2 threshold with full uncertainties. The simulated moments (orange) are shown for comparison.
The covariance matrix of the central moments C 0 is calculated using Gaussian uncertainty propagation C 0 ¼ JCJ ⊺ . Here, J is the Jacobian matrix for the transformation in Eq. (21).
Figure 9 shows the second, third, and fourth central moments as functions of q 2 threshold. The central moments are less correlated with each other than the raw moments but have larger variances. We observe negative correlations between some of the central moments. The full correlation matrix is given in Appendix F. Comparisons of the measured and simulated moments using χ 2 tests show p values greater than 98%.
The Belle Collaboration recently presented a measurement similar to this one [14]. This work provides additional new measurements of the raw and central q 2 moments with comparable precision. We present measurements starting at lower q 2 thresholds of 1.5, 2.0, and 2.5 GeV 2 =c 4 , which retain more information about the inclusive X c spectrum and allow for reductions of the uncertainty on jV cb j. We compare the overlapping measurements of the raw moments from both analyses for q 2 thresholds between 3.0 and 8.5 GeV 2 =c 4 using a χ 2 test including again only measurements with different lower q 2 selections having an observed correlation below 95%. The tests yield p values between 5% and 72%. Here, we assume the systematic uncertainties for the simulation of the X c spectrum are fully correlated between the Belle and Belle II measurements.
We measure the first to fourth moments of the q 2 spectrum of B → X c lν l from 1.5 to 8.5 GeV 2 =c 4 . The precise determinations of these moments are a crucial experimental input for determinations of jV cb j and HQE parameters, proposed by the authors of Ref. [12]. This analysis probes up to 77% of the accessible B → X c lν l phase space, improving on the measurement of Ref. [14], and it includes the experimentally challenging q 2 region of ½1.5; 2.5 GeV 2 =c 4 . The measured moments are also transformed into central moments, which are less correlated but have larger variances than the raw moments.
The uncertainty for the q 2 moments is dominantly systematic, with the uncertainties from the background yield and shape, composition of the X c system, and the simulated detector resolution dominating. A better understanding of the detector and backgrounds will lead to a more precise determination of the q 2 moments in the future and will allow measurements with a q 2 threshold below 1.5 GeV 2 =c 4 .
Recently, a first value of jV cb j was determined using this measurement: Reference [54] finds
which is in good agreement with other inclusive determinations.
We provide numerical results and covariance matrices on HEPData [55]. Spline Smooth Factor 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.01 0.01 0.01 Background Yield and Shape 2.12 1.83 1. 49