Observation of the $B^0_s\!\to D^{*+}D^{*-}$ decay

The first observation of the $B^0_s\!\to D^{*+}D^{*-}$ decay and the measurement of its branching ratio relative to the $B^0\!\to D^{*+}D^{*-}$ decay are presented. The data sample used corresponds to an integrated luminosity of $9\,\text{fb}^{-1}$ of proton-proton collisions recorded by the LHCb experiment at centre-of-mass energies of 7, 8 and $13\,\text{TeV}$ between 2011 and 2018. The decay is observed with more than $10$ standard deviations and the time-integrated ratio of branching fractions is determined to be \begin{align*} \frac{\mathcal{B}(B^0_s\!\to D^{*+}D^{*-})}{\mathcal{B}(B^0\!\to D^{*+}D^{*-})} = 0.269 \pm 0.032 \pm 0.011 \pm 0.008\, , \end{align*} where the first uncertainty is statistical, the second systematic and the third due to the uncertainty of the fragmentation fraction ratio $f_s/f_d$. The $B^0_s\!\to D^{*+}D^{*-}$ branching fraction is calculated to be \begin{align*} \mathcal{B}(B^0_s\!\to D^{*+}D^{*-}) = (2.15 \pm 0.26 \pm 0.09 \pm 0.06 \pm 0.16)\times 10^{-4} \,, \end{align*} where the fourth uncertainty is due to the $B^0\!\to D^{*+}D^{*-}$branching fraction. These results are calculated using the average $B^0_s$ meson lifetime in simulation. Correction factors are reported for scenarios where either a purely heavy or a purely light $B^0_s$ eigenstate is considered.

where the first uncertainty is statistical, the second systematic and the third due to the uncertainty of the fragmentation fraction ratio f s /f d .The B 0 s → D * + D * − branching fraction is calculated to be B(B 0 s → D * + D * − ) = (2.15 ± 0.26 ± 0.09 ± 0.06 ± 0.16) × 10 −4 , where the fourth uncertainty is due to the B 0 → D * + D * − branching fraction.These results are calculated using the average B 0 s meson lifetime in simulation.Correction factors are reported for scenarios where either a purely heavy or a purely light B 0 1 Introduction Decays of B mesons into two charm mesons can be used to probe elements of the Cabibbo-Kobayashi-Maskawa (CKM) matrix [1,2].Measurements of CP violation in B 0 → D ( * )+ D ( * )− , B 0 s → D ( * )+ D ( * )− and B 0 s → D ( * )+ s D ( * )− s decays can be used to determine the CKM angles β [3][4][5][6][7][8] and β s [9], although CP violation in B 0 s → D * + D * − and B 0 s → D * + s D * − s has not yet been measured.These determinations are affected by higherorder Standard Model effects [10][11][12][13][14], which can be constrained using measurements of CP violation parameters or branching ratios in additional decays of the B → DD family, e.g. the B 0 s → D * + D * − decay 1 [15,16].In the B 0 s → D * + D * − decay, tree and penguin transitions cannot contribute.This makes these decays sensitive to the effects of W -exchange and penguin-annihilation diagrams.The dominant Feynman diagrams of the family of B → DD decays are shown in Fig. 1.The absolute branching fraction of the B 0 s → D * + D * − decay is predicted to be (3.1 +1.3 −1.1 ) × 10 −4 using a perturbative QCD approach based on k T factorisation [17].In this paper, the first observation of the B 0 s → D * + D * − decay is presented.Its branching fraction is measured relative to that of the B 0 → D * + D * − decay, thereby canceling systematic uncertainties originating from the uncertainty on the integrated luminosity and the bb cross section.The measurement is performed with data collected by the LHCb experiment, where proton beams collide at centre-of-mass energies up to 1 Charge conjugation is implied throughout the paper.13 TeV.In 2011 and 2012 data samples corresponding to 1 fb −1 and 2 fb −1 were collected at centre-of-mass energies of 7 and 8 TeV, respectively, and from 2015 to 2018 a data sample corresponding to 6 fb −1 was collected at 13 TeV.Due to the different centre-of-mass energies and due to different detector settings, the analysis is performed separately for the two data-taking periods 2011-2012 (Run 1) and 2015-2018 (Run 2), combining both results in the end.
The measurement of the ratio of branching fractions relies on the calculation of the ratio of selection efficiencies and the determination of the signal yields for the two decays.Therefore, the analysis begins with the selection of signal and control mode candidates in Sec. 3 and continues with the extraction of the B 0 s and B 0 yields using mass fits in Sec. 4. The systematic uncertainties are evaluated in Sec. 5 and the final calculation of the ratio of branching fractions is presented in Sec. 6.

Detector and simulation
The LHCb detector [18,19] is a single-arm forward spectrometer covering the pseudorapidity range 2 < η < 5, designed for the study of particles containing b or c quarks.The detector includes a high-precision tracking system consisting of a siliconstrip vertex detector surrounding the proton-proton (pp) interaction region, a large-area silicon-strip detector located upstream of a dipole magnet with a bending power of about 4 Tm, and three stations of silicon-strip detectors and straw drift tubes placed downstream of the magnet.The tracking system provides a measurement of the momentum, p, of charged particles with a relative uncertainty that varies from 0.5% at low momentum to 1.0% at 200 GeV/c.The minimum distance of a track to a primary pp collision vertex (PV), the impact parameter (IP), is measured with a resolution of (15 + 29/p T ) µm, where p T is the component of the momentum transverse to the beam, in GeV/c.Different types of charged hadrons are distinguished using information from two ring-imaging Cherenkov detectors.Photons, electrons and hadrons are identified by a calorimeter system consisting of scintillating-pad and preshower detectors, an electromagnetic and a hadronic calorimeter.Muons are identified by a system composed of alternating layers of iron and multiwire proportional chambers.
Simulation is used to calculate the selection efficiencies and to determine the shapes of the mass distributions.In the simulation, pp collisions are generated using Pythia [20] with a specific LHCb configuration [21].The EvtGen package [22] is used to decay unstable particles, with final-state radiation generated using Photos [23].The interaction of the generated particles with the detector, and its response, are implemented using the Geant4 toolkit [24,25] as described in Ref. [26].

Selection of candidates
The online event selection is performed by a trigger [27,28], which consists of a hardware stage, based on information from the calorimeter and muon systems, followed by a software stage, which applies a full event reconstruction.At the hardware trigger stage, events are required to have a muon with high p T or a hadron, photon or electron with high transverse energy in the calorimeters.The software trigger requires a two-, three-or four-track secondary vertex with a significant displacement from any primary pp interaction vertex.At least one charged particle must have a large transverse momentum and be inconsistent with originating from a PV.A multivariate algorithm [29,30] is used for the identification of secondary vertices consistent with the decay of a b hadron.The B 0 s → D * + D * − and B 0 → D * + D * − decays are reconstructed through D * + → D 0 π + decays, where the D 0 meson decays into the final states K − π + , K + K − or π + π − .These are the dominant two-body D 0 decays with charged final-state particles.Several selection requirements are applied to the final-state tracks, including having a small probability of being a ghost track [31], good track quality and loose particle identification (PID) criteria.The final-state particles and the intermediate D 0 and D * mesons2 are required to have a high χ 2 IP with respect to any PV, where χ 2 IP is defined as the difference in the vertex-fit χ 2 of a given PV reconstructed with and without the candidate being considered.The distance of closest approach between all possible combinations of tracks forming D vertices has to be smaller than 0.5 mm.The cosine of the angle between the momentum of each D meson and the direction from the best PV to the decay vertex has to be greater than zero, where the best PV is the PV for which the B candidate3 has the lowest χ 2 IP .The differences between the reconstructed D * and D 0 invariant masses and their known values [32] have to be smaller than 50 MeV/c 2 and 100 MeV/c 2 , respectively.The reconstructed mass difference between the D * and D 0 candidates is required to be smaller than 150 MeV/c 2 .Both D * mesons must originate from a common vertex, and the B candidate formed must have p T greater than 5000 MeV/c.The B candidates are required to have a χ 2 IP smaller than 25 with respect to the best PV.The decay time of the B meson has to be larger than 0.2 ps to suppress background from particles produced directly in the PV.As a final requirement, a kinematic fit [33] is applied to the decay chain, where the masses of the D mesons and the PV are constrained, and the resulting χ 2 must be less than 50.This fit is used to correctly account for correlations and uncertainties of the vertex positions of the D * and D 0 mesons, particle momenta, flight distances, decay times, and invariant masses.The mass estimate from this constrained fit improves the B mass resolution by around 50%, hence this mass estimate is also used for the extraction of signal candidates in Sec. 4.
A high fraction of events with multiple signal candidates is found in the data, even though the predicted branching fraction is too low for these to be true independent decays.Those are classified into two categories.Multiple candidates that originate from misidentification of a final-state particle occur in 1.4% and 2.5% of all events for Run 1 and Run 2, respectively.These candidates are identified by having the same tracks for all final state particles, but where one K has been misidentified as a π or vice versa.For this category, the candidate corresponding to the decay with the largest branching ratio is chosen.This method is used since it is expected that most of the real decays are present in the D 0 → K − π + channel.Multiple candidates that involve unrelated final-state particles occur in 0.2% and 1.3% of all events for Run 1 and Run 2, respectively, and are rejected randomly so that only one candidate remains per event.decays, to a negligible level.This is validated using simulated samples of these possible backgrounds generated with the RapidSim package [34].

Extraction of signal yields
The total probability density function (PDF) consists of the signal components for the B 0 s → D * + D * − and B 0 → D * + D * − decays and a component for the combinatorial background, i.e., where the PDF and yield of each decay mode is given by P i and N i , respectively.The B 0 s and B 0 signal PDFs both consist of a single double-sided Hypatia function [35].The mean value µ of the B 0 s signal PDF is shifted relative to the mean value of the B 0 PDF by the fixed value µ B 0 s − µ B 0 = 87.26MeV/c 2 calculated with the known meson masses [32].Due to the limited number of selected events in data and simulation, the parameters ζ and β of the generalised hyperbolic distribution G(m − µ, σ, λ, ζ, β) of the Hypatia functions are set to zero to achieve a more stable fit.The simulated samples only contain decays with D 0 → K − π + final states, but the same mass model is used for all final states, as any difference between them is expected to be negligible due to the constraints on both D * masses.The combinatorial background is described with a PDF consisting of an exponential function with a negative slope.The same mass model is used for both data-taking periods with different sets of parameters.
In the fit to the data, the only free parameters are the mean of the B 0 PDF, the widths of both signal PDFs, the slope of the exponential function and the yield of each contribution.The mass distributions and fit projections are shown in Fig. 2 together with the contribution of each component.The resulting yields for the B 0 s → D * + D * − and B 0 → D * + D * − decays are 20 ± 5 and 251 ± 16, respectively, for the data-taking period Run 1.For Run 2 the B 0 s → D * + D * − and B 0 → D * + D * − yields are 79 ± 10 and 1123 ± 34, respectively.For both data-taking periods, the uncertainties are purely statistical.

Systematic uncertainties
The systematic uncertainties are summarised in Table 1.The combined uncertainties are calculated from the individual uncertainties for Run 1 and Run 2, taking into account the correlations.The individual systematic uncertainties are combined in quadrature to obtain the total uncertainties for each data-taking period.
The measurement of the ratio of branching fractions relies on the ratio of b quark hadronisation fractions to B 0 s and B 0 mesons, f s /f d , which is taken from Ref. [36].For the data-taking periods Run 1 and Run 2, all uncertainties of f s /f d are shared between the samples with one additional independent uncertainty of around 0.7% for the Run 1 sample, so overall the uncertainties are partially correlated when the two data-taking periods are combined.Two sources of systematic uncertainty on the efficiency ratio of the B 0 s → D * + D * − and B 0 → D * + D * − modes are considered.The first comes from the limited size of the simulated sample and is treated as uncorrelated between the data-taking periods.Secondly, the nominal efficiency ratio is calculated using simulated unpolarised B 0 s → D * + D * − and B 0 → D * + D * − decays, but these decays involve a pseudoscalar decaying into two vector mesons, which have to be longitudinally or transversely polarised.The effect of different polarisations of the B 0 and B 0 s decay products are analysed by recalculating the efficiency ratio using polarised simulated samples.The polarisation fractions of the B 0 decay are set to their measured values [37].Different polarisation configurations for the B 0 s decay products, ranging from purely longitudinal to purely transverse, are evaluated.The largest deviation to the nominal efficiency ratio is assigned as a systematic uncertainty.The uncertainties of the two data-taking periods are assumed to be fully correlated.
The fraction of events with multiple candidates originating from particle misidentification is of order 2%.In addition to the branching-fraction based method described above, two alternative rejection methods for this kind of multiple candidates are studied to assign a systematic uncertainty: random rejection and rejection according to a probability based on the PID response.For the latter, a ratio for each final state particle of the D 0 decays is calculated indicating how likely the particle is correctly identified.The mass fit is repeated using the alternative rejection methods to extract the new yields.The highest deviation from the default method or between the two alternative methods is assigned as a systematic uncertainty.The uncertainties of the two data-taking periods are assumed to be uncorrelated.
Finally, systematic uncertainties are assigned to account for the choice of the PDF used to describe the signal distributions and to account for the choice of the mass range in the fit.Both uncertainties are estimated using pseudoexperiments.Data sets are generated using alternative mass models of the B 0 and B 0 s contributions, or using alternative mass-fit ranges.The alternative mass model consists of a sum of two Crystal Ball functions [38], which have tails on opposite sides of the peak and parameters that are obtained using fits to simulated samples.The remaining parameters are generated according to their results from the default mass fit.The alternative mass ranges have the upper boundary of the default fit, but the lower boundary is shifted by ±50 MeV/c 2 .The parameters are generated using the default mass fit results for the signal contributions and results of fits in the respective ranges for the background contribution.For each generated data set, an extended maximum likelihood fit is performed using the default model.In the investigation of the alternative mass ranges, the fit uses the respective alternative ranges instead of the default fit range.The deviation of the ratio of the yields relative to the default ratio is calculated for each data set and the mean of the deviations is assigned as the systematic uncertainty.The results are assumed to be fully correlated between the two data-taking periods.

Results
The B 0 s → D * + D * − decay is observed with more than 10 standard deviations, determined using Wilks' theorem [39].
The ratio of branching fractions is given by where ε B 0 and ε B 0 s are the efficiencies of the entire selection, N B 0 and N B 0 s are yields determined from the fit to the D * + D * − mass distribution, and f s /f d is the ratio of hadronisation fractions of B 0 and B 0 s , which is taken from Ref. [36].Accordingly, the f s /f d ratios are given by Calculating the selection efficiencies using simulated samples results in The ratios of the numbers of signal candidates resulting from the fit are calculated to be where the first uncertainty is statistical, the second is systematic and the third is due to the uncertainty of the fragmentation fraction ratio f s /f d .The results are combined by calculating a weighted average, using the total uncertainties on the individual measurements and taking into account correlations.The combined ratio of branching fractions is Using the known value of the branching fraction of the B 0 → D * + D * − decay [32], the B 0 s → D * + D * − branching fraction is calculated to be where the fourth uncertainty is due to the B 0 → D * + D * − branching fraction.The central value of the time-integrated branching fraction presented in this paper depends on the lifetime of the B 0 s meson [40].There is currently no measurement or prediction for the B 0 s lifetime in this decay, so the average B 0 s lifetime is used in the generation of simulated B 0 s → D * + D * − decays.However, the light and heavy eigenstates have significantly different lifetimes, τ (B 0 s,L ) = 1.423 ps and τ (B 0 s,H ) = 1.620 ps [32].Since the B 0 s selection efficiency depends on the lifetime, correction factors for the efficiency are calculated for the scenarios where either a purely heavy or a purely light B 0 s eigenstate is considered.The time-integrated correction factors for the efficiency are 0.947 and 1.045 for light and heavy B 0 s eigenstates, respectively.If the lifetime of the B 0 s meson in the B 0 s → D * + D * − decay is found in further studies to differ significantly from the average B 0 s lifetime, the central value of the paper result can be corrected.The equivalent effect also occurs for the B 0 → D * + D * − decay but is neglected due to the small value of ∆Γ d [32].

Summary
The B 0 s → D * + D * − decay is observed with more than 10 standard deviations.Its branching fraction is measured relative to the B 0 → D * + D * − branching fraction.The measurement is performed using data collected by the LHCb experiment in Run 1 and Run 2 corresponding to an integrated luminosity of 9 fb −1 .The time-integrated ratio of branching fractions is measured to be where the first uncertainty is statistical, the second systematic and the third is due to the uncertainty of the fragmentation fraction ratio f s /f d .The ratio is measured assuming an average B 0 s lifetime in the calculation of the B 0 s selection efficiency, since no measurement or prediction for the B 0 s lifetime is available in this decay.Correction factors for the efficiency are reported in Sec.6 for the scenarios where either a purely heavy or a purely light B 0 s eigenstate is considered, so that the central value of the result presented here can be recalculated when a measurement or prediction is available.The B 0 s → D * + D * − branching fraction is determined to be B(B 0 s → D * + D * − ) = (2.15 ± 0.26 ± 0.09 ± 0.06 ± 0.16) × 10 −4 , where the fourth uncertainty is due to the B 0 → D * + D * − branching fraction.The result is in agreement with the theoretical prediction using a perturbative QCD approach based on k T factorisation [17].It can be used to constrain higher-order SM effects in CP violation observables in B 0 → D ( * )+ D ( * )− and B 0 s → D

Figure 1 :
Figure 1: Dominant Feynman diagrams for B → DD decays comprising (top left) tree, (top right) penguin, (bottom left) W -exchange and (bottom right) penguin-annihilation transitions.The B 0 s → D * + D * − decay can only occur via the W -exchange or penguin-annihilation diagram.

Figure 2 :
Figure 2: Distributions of the D * + D * − mass for (left) Run 1 and (right) Run 2 data samples, with the fit results overlaid.In addition to the full fit result (solid black), each of the individual components are shown: the B 0 → D * + D * − signal (dotted blue), the B 0 s → D * + D * − signal (dash-dotted red) and the combinatorial background (dashed green).

Table 1 :
Systematic uncertainties relative to the measured value of the ratio of branching fractions for Run 1, Run 2, and their combination, taking into account correlations.