Prospects for $B_{c}^+\to \tau^+ \nu_\tau$ at FCC-ee

This paper presents the prospects for a precise measurement of the branching fraction of the leptonic $B_{c}^+\to \tau^+ \nu_\tau$ decay at the Future Circular Collider (FCC-ee) running at the $Z$-pole. A detailed description of the simulation and analysis framework is provided. To select signal candidates, two Boosted Decision Tree algorithms are employed and optimised. The first stage suppresses inclusive $b\bar{b}$, $c\bar{c}$, and $q\bar{q}$ backgrounds using event-based topological information. A second stage utilises the properties of the hadronic $\tau^{+} \to \pi^+ \pi^+ \pi^- \bar{\nu}_\tau$ decay to further suppress these backgrounds, and is also found to achieve high rejection for the $B^+\to \tau^+ \nu_\tau$ background. The number of $B_{c}^+\to \tau^+ \nu_\tau$ candidates is estimated for various Tera-$Z$ scenarios, and the potential precision of signal yield and branching fraction measurements evaluated. The phenomenological impact of such measurements on various New Physics scenarios is also explored.


Introduction
Leptonic pseudoscalar meson decays such as B + c →τ + ν τ are theoretically clean probes to test for the presence of physics beyond the Standard Model (SM). The only hadronic inputs required to compute their decay branching fractions in the SM are the decay constants, which have been precisely determined for several transitions by means of numerical simulations of QCD on the lattice (LQCD) [1]. In the past several years, numerous discrepancies from SM predictions have been observed in tree-level [2][3][4][5][6][7][8][9] and loop-induced [10][11][12][13] semileptonic b-hadron decays, often referred to as the B-physics anomalies. The B + c →τ + ν τ decay 1 can be directly related to the anomalies in tree-level decays as they occur through the same quark-level transition, b → cτ ν τ , thus offering a clean and independent check of these experimental results [14,15]. Furthermore, B + c →τ + ν τ decays are highly sensitive probes of pseudoscalar contributions from New Physics (NP), as predicted for instance in extensions of the SM Higgs sector, such as Two-Higgs-Doublet Models (2HDM) [16], as well as in specific leptoquark models [17,18].
Despite the fairly large branching ratio of the B + c → τ + ν τ decay within the SM (≈ 2%), an observation and precise measurement of its properties are very challenging at a hadron collider. This is due to the large missing energy in the final state, no knowledge of the centre-of-mass energy of the bb production process, high backgrounds due to the hadronic environment with multiple primary vertices (PVs), and the lack of any reconstructible B + c decay vertex. Moreover, these decays cannot be studied at B-factories since B + c mesons are too heavy to be produced from Υ(5S) → bb decays. For these reasons, future Z-factories provide a unique environment to study these processes in the future. Indeed, the large number Z bosons produced, up to N Z ≈ 5 × 10 12 in the case of FCC-ee (so-called Tera-Z), together with the possibility to constrain the missing energy from the neutrinos, would make this measurement possible. In this paper, the feasibility of performing such a measurement at FCC-ee in Z-pole operation [19][20][21] is demonstrated, employing the τ + → π + π + π −ν τ mode to provide the τ + decay vertex and thus a measure of the combined B + c and τ + flight distance. This reconstruction method offers additional means of background rejection, due to the lower lifetime of the B + c meson relative to lighter b-hadrons, and the different resonant properties of hadronic τ + decays compared to backgrounds from b-and c-hadrons. In the context of the CEPC project, a similar study has been presented in Ref. [22], where leptonic decays of the τ + were considered.
The branching fraction of the B + c →τ + ν τ decay in the SM can be written as where G F is the Fermi constant, m Bc and m τ are the the B + c meson and τ + lepton masses, respectively, τ Bc denotes the B + c meson lifetime [22], V cb is the CKM matrix element for b → c transitions [23], and f Bc represents the B + c meson decay constant, f Bc = 427 (6) MeV, which has been computed via LQCD simulations in [24] (see also Ref. [25]). By combining these inputs with the latest value of |V cb | excl. = 39.09(68) × 10 −3 [1], determined from exclusive B → D ( * ) ν decays using the BGL parameterization for the B → D ( * ) form factors [26,27], is obtained. One of the challenges of a precise measurement of the B + c →τ + ν τ decay branching fraction is to properly normalise the measurement, in order to avoid relying on the unknown B + c meson hadronisation fraction, f (B ± c ) ≡ f (b → B ± c ) [28]. This was the main caveat of previous attempts to extract limits on B(B + c → τ + ν τ ) from LEP data [29,30], as discussed for example in Refs. [31,32]. In this work, the possibility of normalising the measurement to the B + c → J/ψµ + ν µ mode is proposed, using the B + c → J/ψ form factors recently computed via LQCD in Ref. [33,34] to predict B(B + c → J/ψµ + ν µ ). The B + c → τ + ν τ branching fraction is then determined as follows, where the number of signal and normalisation candidates, N (B + c → τ + ν τ ) and N (B + c → J/ψµ + ν µ ), are both measured in data.
Their corresponding total efficiencies (B + c → τ + ν τ ) and (B + c → J/ψµ + ν µ ) are estimated using numerical simulations. The branching fractions of the J/ψ to a pair of muons and the τ + hadronic decay to three pions are taken from Ref. [35], and the B + c → J/ψµ + ν µ branching fraction can be accurately predicted by using the LQCD results for the form-factors to be B(B + c → J/ψµ + ν µ ) SM = 0.0135 (11) [33,34], where the |V cb | excl. value quoted above is used. Alternatively one could use the fully reconstructed B + c → J/ψπ + decay mode as a normalisation instead of the partially reconstructed B + c → J/ψµ + ν µ decay. However, the branching fraction of the B + c → J/ψπ + decay is measured to be around 5% of the B + c → J/ψµ + ν µ branching fraction [36], and as such the uncertainty due to finite normalisation statistics would be considerably larger. In addition, no reliable SM prediction for B(B + c → J/ψπ + ) is yet available. Taking B + c → J/ψµ + ν µ as the normalisation mode, the ratio can be measured. This ratio also offers a sensitive probe of NP contributions, with the advantage of being independent of both the B + c production rate and |V cb |, as discussed in Refs. [37][38][39] for example. The limiting factor on the precision of R c is the uncertainty on the B + c → J/ψ form-factors, which amount to a theory uncertainty of ≈ 7% on R c . This source of uncertainty can be improved in the future with updated LQCD computations, the precision of which can be considerably improved by performing a dedicated angular analysis of B + c → (J/ψ → µ + µ − )µ + ν µ decays with the upgraded LHCb detector or at FCC-ee. Such measurements would be particularly useful to constrain the form factors in the large-recoil region.
The remainder of this paper is organised as follows: a description of the experimental setup, software, and simulated samples is provided in Section 2; a demonstration of the multivariate selections used to select signal decays with high purity is provided in Section 3, along with estimates of the achievable signal yield and branching fraction precision as a function of N Z ; the phenomenological impact of this measurement on well-motivated NP scenarios is discussed in Section 4.
The work in this paper has been conducted after the publication of the FCC conceptual design reports [20,21,40] and after the European Strategy Update for Particle Physics released its recommendation to investigate the technical and financial feasibility of a future hadron collider at CERN with a centre-of-mass energy of at least 100 TeV and with an electron-positron Higgs and electroweak factory as a possible first stage. This article provides a detailed description of the key ingredients for physics analyses at FCC-ee.

FCC-ee
The international Future Circular Collider (FCC) study aims at a design of p-p, e + e − , and e-p colliders to be built in a new 100 km tunnel in the Geneva region. The e + e − collider (FCC-ee) has a centre of mass energy range between 91 (Z-pole) and 375 GeV (tt). The FCC-ee offers unprecedented possibilities for measuring the properties of the four heaviest particles of the SM (the Higgs, Z, and W bosons, and the top quark), but also those of the b and c quarks and of the τ lepton. In addition, circular colliders have the advantage of delivering collisions to multiple interaction regions, which allow different detector designs to be studied and optimised -up to four are under consideration for FCC-ee. Moreover, the huge statistics anticipated at the Z peak (the so-called "Tera-Z" run) brings specific challenges, as the systematic uncertainties of the measurements should be commensurate with their small statistical uncertainties.

Simulation of the detector response
The detector response has been simulated via the DELPHES software package [41]. It is a C++ framework, performing a fast multipurpose detector response simulation. The simulation includes a tracking system embedded in a magnetic field, calorimeters, and a muon system. The framework is interfaced to standard file formats (e.g. Les Houches Event File or HepMC) and outputs observables such reconstructed charged tracks which can be used for dedicated analyses. The simulation of the detector response takes into account the effect of the magnetic field, the granularity of the calorimeters, and sub-detector resolutions. In the pre-release 3.4.3pre10 used for this analysis, DELPHES provides parameterised track information with the full covariance matrix using the FastTrackCovariance software.
The detector configuration considered is the Innovative Detector for Electron-positron Accelerators (IDEA) concept. It comprises a silicon pixel vertex detector, a large-volume extremely-light short-drift wire chamber surrounded by a layer of silicon micro-strip detectors, a thin, low-mass superconducting solenoid coil, a pre-shower detector, a dual-readout calorimeter, and muon chambers within the magnet return yoke [21]. The DELPHES configuration card used for this analysis is accessible in the repository given in Ref. [42]. Finally, the k4SimDelphes [43] project converts DELPHES objects to EDM4hep [44], and the subsequent Monte Carlo (MC) production is performed in the common EDM4hep data format.

Monte-Carlo Production
MC event samples are used to simulate the response of the detector to signal and background processes. Signal and background events are generated with Pythia [45] version 8.303 using the leading order cross-section from the generator with no K-factor. Decays of unstable particles are described using EvtGen [46] version 02.00.00, in which final-state radiation is generated using Photos [47]. Several important parameters are configured to be the same for all samples in the Pythia steering cards. At FCC-ee [21], the energy of the beams is distributed according to a Gaussian function. At the Z peak, the beam energy spread amounts to 0.132% of the incoming beam energy, half the Z boson mass, which equates to 0.0602 GeV. The position of the interaction region depends on the running conditions of the machine. At the Z-pole the bunch length is σ z = 12.1 mm. The bunch dimensions in the transverse plane, at the interaction point (IP), are given by σ x,y = β * x,y × x,y , where the values of the β function at the IP, and the horizontal and vertical emittance x,y are given in the FCC-ee Conceptual Design Report [21] and are approximately σ x = 6.4 µm and σ y = 28.3 nm. For Gaussian bunches, the PV distribution in (x, y, z) is well approximated by a 3-dimensional Gaussian distribution, with where α denotes the half-crossing angle of 15 mrad. This yields values for the primary vertex smearing, σ P V x = 4.5 µm, σ P V y = 20 nm, and σ P V z = 0.3 mm. All of the generator configurations and steering cards are documented and preserved [42]. The production of Monte-Carlo events is achieved using the FCC common tools and CERN computing and storage resources [48]. Approximately 10 10 events are produced, representing about 55 TB of disk space. Dedicated productions with orthogonal seeds between the analysis and multivariate training samples have been considered in order to avoid over-training.

Analysis framework
A sophisticated analysis framework has been developed for all FCC analyses using the common EDM4hep data format. It is based on RDataFrames [49], where C++ code is conveniently compiled in a ROOT [50] dictionary as "analysers" which are subsequently called in Python [51]. Several external packages such as ACTS [52], FastJet [53], and awkward [54] are included. The analysis code is distributed via the CERN virtual machine file system cvmfs, and can be run locally or on batch systems. The complete software stack used to produce the results in this paper can be accessed and the results reproduced [55].

Specificity for this analysis
For this analysis, dedicated new features of the analysis framework have been developed and are explained in this section.
Perfect vertex seeding: Excellent vertex finding is crucial for this analysis. While detailed investigations are ongoing to estimate the impact of imperfect vertexing, for the following results it is assumed that vertices can be perfectly seeded. The procedure is to first find all of the MC vertices by selecting stable charged particles originating from the same point. From those MC vertices, reconstruction-level vertices are fitted using the reconstructed tracks associated to the MC particles attached to the MC vertex. A plot comparing the MC and reconstruction-level number of vertices is shown in Fig. 1 (a). This procedure properly takes into account the migration of higher number of MC tracks to a given reconstructed multiplicity, as illustrated in Fig. 1 (b), where about 7% of the three-track reconstructed vertices originate from a four-track MC vertex.
Perfect particle identification: In the energy range considered for the identification of pions and kaons in the analysis (∼ 10 GeV), extremely good discrimination is expected, thus for this measurement it is considered that the pions and kaons can be perfectly identified.
Processing the first stage of the analysis over the full sample statistics, and calculating complex quantities (such as the thrust described in Section 3.2, vertexing, building the candidates), takes approximately half a day on a batch system. The total sample size after first-stage processing is ∼ 280 GB, representing a reduction factor of 200. The second stage of the analysis can then be run locally very quickly, as can all of the final analysis steps described in Section 3.

Analysis
To demonstrate the feasibility of a measurement of the B + c → τ + ν τ mode at FCC-ee, a selection procedure based on the differences in reconstructed event and candidate properties in simulated B + c → τ + ν τ and inclusive hadronic Z decays is developed. The selection consists of a series of rectangular cuts in addition to two boosted decision tree (BDT) classifiers, which together achieve a high purity final dataset with a clearly identifiable signal component. The selection is now described, and estimates for the signal yield precision derived via template fits to combined samples of signal and background decays. A discussion of the potential branching ratio precision as a function of N Z is also given. All the tools used for the following analysis are available in Ref. [56].

Signal and background samples
All samples used in the selection studies are generated and processed using the framework detailed in Sec. 2. For the BDT training, orthogonal samples of signal and background events are generated; these samples are not used in any other analysis steps to avoid biasing the BDT distributions. Simulated samples of B + c → τ + ν τ with τ + → π + π + π −ν τ decays are considered as signal in the selection studies. The samples are generated using EvtGen, where the B + c is decayed using the SLN model and the τ + → π + π + π −ν τ decay is generated using the TAUHADNU model. The TAUHADNU model is used instead of alternative models such as TAUOLA to enable highly efficient truth-matching of the pions to a τ + parent; the difference in τ + decay product kinematics across models is not sufficient to alter the outcome of the selection studies.
For both BDT training stages, inclusive samples of Z → bb, cc, and qq decays are used as background, where q ∈ {u, d, s}. The samples are generated using Pythia, and are found to have consistent distributions when compared to inclusive samples generated using EvtGen. The background samples are combined according to known hadronic Z branching fractions [35] and the total efficiencies of selection cuts applied prior to the training steps.
Prior to the BDT cut optimisation, none of the 10 9 inclusive Z → cc and Z → qq events are found to pass sufficiently tight cuts on both BDTs. As such, background from these sources is not considered in the optimisation or subsequent fit studies. After the same cuts, the remaining statistics in the inclusive Z → bb sample are found to be insufficient for determining the background rejection accurately in the cut optimisation. To boost the background statistics for the optimisation, samples of exclusive b-hadron decays are generated, where the decay modes are chosen based on the composition of the remaining inclusive Z → bb sample. The following decays are considered: In each of the exclusive b-hadron samples, all of the b-hadron decay products are decayed inclusively. The list of exclusive decays considered is not exhaustive, and covers around 10% of the decay width for each B hadron. As a result, a factor 2.5 difference in rate relative to the inclusive Z → bb sample is observed after tight BDT cuts. This factor is used to scale the exclusive sample yield estimates in the optimisation procedure, in order to avoid underestimating the expected background level.

Thrust axis and event hemisphere definitions
The signal selection relies on the large missing energy signature of B + c → τ + ν τ decays, which arises due to the presence of two neutrinos in the final state. In a Z → bb event involving a signal decay, the signal side of the event will on average contain considerably more missing energy than the non-signal side. This is in contrast to a general Z → bb event, where both sides of the event will contain similar amounts of reconstructed energy on average. To determine the energy imbalance in an event, it is necessary to divide the event into two hemispheres, each corresponding to one of b-quarks produced in the Z decay.
In this analysis, hemispheres are defined event-by-event using the plane normal to the thrust axis. The thrust is the unit vectorn which minimises where p i is the momentum vector of the i th reconstructed particle. The axis along whichn lies is referred to as the thrust axis, and provides a measure of the direction of the quark pair produced in the Z decay. Reconstructed particles are assigned to either hemisphere based on the angle θ between their momentum vector and the thrust axis. In this analysis, the thrust is defined to point towards the hemisphere with less total energy; for a particle in the minimum energy hemisphere, cos(θ) ≥ 0, while particles in the maximum energy hemisphere have cos(θ) < 0.

First-stage BDT
The first step of the selection is designed to separate signal and background decays based on the energy signatures of both hemispheres and other general properties of the event. Prior to the first-stage BDT, events are required to contain a reconstructed PV and at least one reconstructed 3π candidate with an associated vertex. In addition, at least one of the 3π candidates is required to reside in the minimum energy hemisphere, since this hemisphere is more likely to contain the signal decay.
The first-stage BDT is trained using xgboost [57] with a sample of 7 × 10 5 signal events passing the above pre-selection and a combined sample of one million inclusive Z → bb, cc, and qq events. The Z → bb background sample is filtered to remove B + c → τ + ν τ and B + → τ + ν τ events. The relative proportion of each Z decay type in the background sample is determined using their known hadronic Z branching fractions [35], multiplied by pre-selection efficiencies determined by the ratio of the number of events passing the pre-selection cuts relative to the number of generated events. The BDT is trained using the following features: • Total reconstructed energy in each hemisphere; • Total charged and neutral reconstructed energies in each hemisphere; • Charged and neutral particle multiplicities in each hemisphere; • Number of tracks in the reconstructed PV; • Number of reconstructed 3π candidates in the event; • Number of reconstructed vertices in each hemisphere; • Minimum, maximum, and average radial distance of all decay vertices from the PV.
The performance of the BDT is illustrated in Fig. 2, where the BDT distributions in signal, B + → τ + ν τ , and each category of inclusive Z background are shown alongside their corresponding efficiency profiles. The BDT is found to have a ROC curve area of 0.984, highlighting the excellent rejection of inclusive background achieved. Although B + → τ + ν τ decays are not considered in the training background, the BDT achieves some rejection of this mode relative to signal. This is due to the different event-level properties of B + → τ + ν τ and B + c → τ + ν τ decays, which arise since the B + c meson is produced with an associated charm quark that results in production of an associated charm hadron. Due to the finite charm hadron lifetime, this results in less energy and fewer tracks at the primary vertex on average compared to a B + → τ + ν τ event, as well as more reconstructed displaced vertices.

Second-stage BDT
The second stage of the selection focuses on properties of the reconstructed 3π candidate and properties of other reconstructed decay vertices in the event. Prior to training the second-stage BDT, events are required to pass a cut of > 0.6 on the first stage BDT. This cut is over 90% efficient on signal, and removes more than 90% of all background types. As shown in Fig. 2, the Z → bb background is rejected least by the first-stage BDT, since it predominantly involves both b → cW and c → sW quark transitions, leading to more missing energy in the case of leptonic W decays. In addition to the first-stage BDT cut, the difference in energy between the maximum and minimum energy hemispheres is required to exceed 10 GeV/c 2 , in order to retain more signal-like events with a large energy imbalance.
In each event, a single 3π candidate is chosen as the signal candidate. The signal candidate must reside in the minimum energy hemisphere, and must have the smallest vertex fit χ 2 of all 3π candidates in that hemisphere. Selected 3π candidates are required to have an invariant mass below that of the τ lepton, and must have at least one m(π + π − ) combination within the range 0.6 − 1.0 GeV/c 2 . These cuts retains candidates consistent with the a 1 (1260) + → (ρ 0 → π + π − )π + decay, via which all τ + → π + π + π −ν τ decays proceed. The second-stage BDT is also trained using xgboost, with a sample of 5 × 10 5 signal events and a combined sample of one million inclusive Z decays. The proportions of Z → bb, cc, and qq decays in the combined sample are determined using their known Z branching fractions and the efficiencies of all pre-selection cuts.
The BDT is trained on the following features: • 3π candidate mass, and masses of the two π + π − combinations; • Number of 3π candidates in the event; • Radial distance of the 3π candidate from the PV; • Vertex χ 2 of the 3π candidate; • Momentum magnitude, momentum components, and impact parameter (transverse and longitudinal) of the 3π candidate; • Angle between the 3π candidate and the thrust axis; • Minimum, maximum, and average impact parameter (longitudinal and transverse) of all other reconstructed decay vertices in the event; • Mass of the PV; • Nominal B energy, defined as the Z mass minus all reconstructed energy apart from the 3π candidate.
The performance of the BDT is illustrated in Fig. 3, where the BDT distributions in signal, B + → τ + ν τ , and each category of inclusive Z background are shown alongside their corresponding efficiency profiles. The BDT is found to have a ROC curve area of 0.966, indicating the high rejection of background achieved even after the BDT1 > 0.6 cut. The BDT is also found to reject B + → τ + ν τ decays to a high level, owing to the larger lifetime of the B + meson compared to the B + c meson which results in a greater 3π displacement from the PV on average. In addition, the lack of an associated charm hadron in B + → τ + ν τ decays is also discriminated against by the second-stage BDT. To select a high-purity sample of signal decays, an optimisation procedure is employed to tune the two BDT cuts, using estimates for the signal and background yields expected at FCC-ee.

Selection optimisation
To determine the best BDT cut values, a two-dimensional optimisation procedure is performed. Prior to the optimisation, cuts of BDT1 > 0.99 and BDT2 > 0.99 are applied in order to focus on the signal region. A grid of 50 cuts for both BDT1 and BDT2 between 0.99 and 1.0 is scanned (2500 points in total), and at each point the expected signal (S) and background (B) yields are estimated. The point where the signal purity P = S/(S + B) is maximised is taken to represent the best cut values for BDT1 and BDT2.
To estimate the signal yield at a given pair of BDT cuts, the following formula is used: (6) where N Z is the number of Z bosons produced at FCC-ee, B(Z → bb) is the known branching fraction of the Z → bb decay [35], the factor of two accounts for the fact that either b-quark from the Z decay can produce a B + c meson, f (B + c ) = 0.04% is the B + c meson production fraction taken from Pythia, B(B + c → τ + ν τ ) = 1.94% is the SM prediction for the signal decay branching fraction, B(τ + → π + π + π −ν τ ) is the known τ + → π + π + π −ν τ branching fraction [35], and is the efficiency of the full selection at a given pair of BDT cuts determined from simulation.
With sufficiently tight cuts to BDT1 and BDT2, all 10 9 generated inclusive Z → cc and Z → qq decays are rejected. As such, background from these sources is not considered in the optimisation routine. Insufficient statistics remain at tight BDT cuts in the inclusive Z → bb sample to measure accurate background rejection figures. As such, the exclusive background samples described in Sec. 3.1 are used to represent the remaining sources of background. The yield for a particular exclusive background decay B → DX at a particular pair of BDT cuts is given by where f (B) is the hadron production fraction for a particular b-hadron taken from Pythia (f is the decay mode branching fraction taken from Ref. [35] where measured, and is the efficiency determined using simulation. Where background mode branching fractions are not yet measured, assumptions are made based on the measured branching fractions of the most topologically similar decay modes. The total background level is given by a sum over all exclusive modes considered, with a multiplicative factor of 2.5 to account for the observed difference in rate between the inclusive Z → bb and exclusive b-hadron samples. To determine the background efficiencies, per-mode efficiencies for a combined BDT1 > 0.95 and BDT2 > 0.95 cut are first determined, as sufficient statistics remain in each exclusive sample to measure these efficiencies accurately. The efficiencies for subsequent BDT1 and BDT2 cuts relative to this point are then measured using a combined sample of all exclusive decay modes. The distributions above 0.95 are parameterised using cubic spline functions s 1 (x 1 ) and s 2 (x 2 ), where x 1 and x 2 represent the BDT1 and BDT2 values. The combined efficiency for x 1 > α and x 2 > β cuts is given by where m 1 and m 2 are the maximum BDT1 and BDT2 scores observed in the summed background sample, respectively. The total efficiency for each background mode is then given by = ( The spline fits to the summed exclusive background BDT distributions are shown in Fig. 4, where an example cut of > 0.99 is shown along with the optimal cuts of BDT1 > 0.99979 and BDT2 > 0.99693 found by the optimisation procedure. At these optimal cut values, the following yields are estimated for N Z = 5 × 10 12 : where the expected B + → τ + ν τ yield is calculated using the measured branching fraction for this mode [35]. The total signal efficiency is found to be 0.39%, and the signal purity is determined to be 85%, demonstrating the excellent performance of the two BDTs in reducing the background from b-hadron decays. The 6.6% rate of the B + → τ + ν τ mode relative to signal is also notable; given the factor 10 3 higher production rate of B + mesons relative to B + c mesons according to Pythia, and a relative branching ratio B(B + → τ + ν τ )/B(B + c → τ + ν τ ) ∼ 0.5%, the B + mode is expected to contribute at five times the level of the signal decay prior to any selection cuts. The full selection efficiency for the B + → τ + ν τ mode is found to be 4.3 × 10 −5 , owing to the excellent rejection of BDT2 in particular. The efficiencies for each exclusive background mode considered are given in App. A Table 2; efficiencies at the 10 −10 − 10 −9 level are found.

Fit to measure the signal yield
To evaluate the potential precision of a signal yield measurement with N Z = 5 × 10 12 , pseudoexperiment fit studies are performed. To select a fit variable, comparisons between all signal and background variable distributions are performed after tight BDT cuts; only those variables related to hemisphere energy are found to provide discrimination. Of all considered variables, the total energy in the maximum energy hemisphere is found to provide the most discrimination. Normalised distributions of this variable in signal events, B + → τ + ν τ events, and exclusive background events are shown in Fig. 5 (a), where cuts of > 0.99 are applied to both BDTs. The exclusive background distribution shown is a sum of all exclusive b-hadron modes considered, where the modes are combined according to their expected yields from the cut optimisation.
In the selection, the signal 3π candidate is required to reside in the minimum energy hemisphere, since this is most likely to be the true signal hemisphere due to the large missing energy. As a result, in signal events the energy distribution in the maximum energy hemisphere closely resembles an inclusive b-quark decay from a Z boson; a peak close to m(Z)/2 with a tail extending to lower energies due to missing energy. A very similar distribution is also observed in B + → τ + ν τ events. The total charged and neutral energies in the maximum energy hemisphere are found to correspond closely in B + (c) → τ + ν τ events, as shown in Fig. 5 (b). The narrow diagonal observed corresponds to the peaking structure in total hemisphere energy close to m(Z)/2, as shown in Fig. 5 (a).
In Z → bb background events, the distinction between the minimum and maximum energy hemispheres is not well defined, as either side of the event could contain more missing energy. Prior to any selection cuts, both hemispheres exhibit the inclusive structure with a peak close to m(Z)/2 and a tail. However, after application of the full selection, the energy distribution in the maximum energy hemisphere is biased downwards and the peaking structure is no longer observed. Both the charged and neutral energies in the maximum energy hemisphere are biased downwards by the selection, as shown in Fig. 5 (c). The weaker correspondence between the neutral and charged energies along the diagonal results in no clear peaking structure and a spread to lower values, as shown in Fig. 5 (a). These changes to the energy distributions are found both in the inclusive Z → bb sample and in the summed exclusive b-hadron sample. For the pseudoexperiment fit studies, the exclusive sample is used in order to retain sufficient statistics for the creation of a template.
To create a total probability density function (PDF) for the pseudoexperiment fit, the normalised histogram templates shown in Fig. 5 are each multiplied by the corresponding yield estimates from the cut optimisation and then summed. To create a pseudoexperiment dataset, a copy of the total PDF is created and each bin varied independently according to Poisson statistics. A fit is then performed to the pseudoexperiment dataset using the total PDF, where the yield of the signal and background components are free parameters. As the B + → τ + ν τ distribution is very similar to signal, and this mode contributes at only the ∼ 7% level relative to signal, it is not possible to freely vary N (B + → τ + ν τ ) in the fit. As such, the yield of this component is constrained according to a Gaussian, with a central value of 249 events (the expected yield from the optimisation) and a width of 10 events. This width corresponds to 5% relative uncertainty on the B + → τ + ν τ yield, which is the anticipated precision on B(B + → τ + ν τ ) from Belle II [58]. Given the small anticipated contribution from the B + → τ + ν τ mode relative to signal, the absolute uncertainty on B(B + → τ + ν τ ) will not contribute a significant source of uncertainty to N (B + c → τ + ν τ ). The result of a single pseudoexperiment fit is shown in Fig. 6 (a), and the distribution of signal yields measured in 2000 fits to different generated pseudoexperiment datasets is shown in Fig. 6 (b). Both the signal and background yields are measured without bias in either their central values or uncertainties. From a fit to the distribution of signal yields, the signal yield uncertainty is found to be 101, which corresponds to a relative uncertainty of 2.4% on the generated yield of 4295 events. This uncertainty is purely statistical, and does not account for potential sources of systematic uncertainty. However, the expected 5% relative uncertainty on B(B + → τ + ν τ ) is included via the Gaussian constraint on N (B + → τ + ν τ ) in the fits. Fits performed with ten times more background included in the toy samples are also found to be stable, with a relative signal yield uncertainty of 2.9%. The analysis is thus robust to large increases in background yield relative to what is modelled here.

Fit performance for different N Z
The fit results shown in Fig. 6 correspond to a final sample selected from a dataset containing N Z = 5×10 12 . The fit performance is also studied at lower sample sizes of N Z = [0.5, 1, 2, 3, 4]×10 12 , in order to evaluate the signal yield precision possible at earlier stages of FCC-ee operation. For each N Z value, the cut optimisation is rerun in order to maximise the purity; highly consistent optimal purity is found across all N Z values. Sets of 2000 pseudoexperiment fits are run for each N Z value, using the expected signal, B + → τ + ν τ , and background yields from the cut optimisation. The signal yields expected as a function of N Z , as well as their uncertainties as measured in the pseudoexperiment fits, are summarised in Tab. 1. The relative signal yield precision as a function of N Z is illustrated in Fig. 7, where four different systematic uncertainty scenarios are shown; σ syst = [0, 0.25, 0.5, 1] × σ stat . All values shown are summarised in App. B Tab. 3. The level of systematic uncertainty in a real analysis will depend on several factors, such as: • Detector resolution, reconstruction efficiency, and calibration quality; • The size of the simulated samples used to create fit templates; • The decay models used to generate signal and background decays; • Knowledge of the relative proportions of decay modes entering the total background template.
Given the high signal purity achievable, however, and the distinctive shape of the signal maximum hemisphere energy distribution, an eventual measurement is not expected to be limited by systematic uncertainties. Assuming that the systematic uncertainties can be controlled at the level σ syst = σ stat , the relative precision possible on N (B + c → τ + ν τ ) with N Z = 5 × 10 12 is

Branching fraction determination
It is common to measure signal modes relative to a normalisation decay, in order to minimise systematic uncertainties and cancel the effects of hadron production. One suitable choice of normalisation mode for B + c → τ + ν τ is the semileptonic B + c → J/ψµ + ν µ decay, where the J/ψ → µ + µ − channel can be selected in order to provide a clean three-muon B + c decay vertex. This mode can be reconstructed and selected with high efficiency, as sources of lighter b-hadron background can be eliminated with a m(J/ψµ) > 5.3 GeV/c 2 cut [36,59]. Above this cut, the only sources of remaining   background are from random combinations of three muons (expected to be small at a FCC-ee) and contributions from B + c → J/ψµ + ν µ X decays where X is not considered in the invariant mass sum. The latter contribution, from decays such as B + c → (ψ(2S) → J/ψπ + π − )µ + ν µ , can be reduced using isolation requirement, where all other charged particles and neutrals in the signal hemisphere must be inconsistent with originating from the 3µ vertex.
With B + c → J/ψµ + ν µ as a normalisation mode, the ratio of branching fractions can be measured. Assuming SM amplitudes in the normalisation decay, the ratio R c is highly sensi-tive to NP couplings to τ leptons, and is independent of the value of |V cb | and the B + c hadronisation fraction f (B + c ). The J/ψ and τ branching fractions are well measured [35], while the signal and normalisation efficiencies can be determined with high accuracy given sufficiently large simulated samples. As the signal and normalisation modes both involve the reconstruction and selection of three charged tracks from a common vertex, systematic uncertainties in the absolute efficiencies are expected to cancel to a high degree in the ratio. The SM prediction for B(B + c → J/ψµ + ν µ ) SM = 0.0135 ± 0.0011, which is calculated using lattice QCD and the current value of |V cb | excl. [1]. Assuming a selection efficiency of 10% for this mode, the anticipated normalisation yield is around 50, 000 events with N Z = 5 × 10 12 . The normalisation yield is thus an order of magnitude larger than the signal yield, and as such will not contribute significantly to the branching fraction ratio uncertainty. Thus, the precision on R c is expected to be dominated by the uncertainty on N (B + c → τ + ν τ ). The anticipated relative uncertainty on R c as a function of N Z is shown in Fig. 8 (a), where the signal yields and uncertainties from Tab. 1 are used as input. The uncertainties shown also include the current uncertainties on the J/ψ and τ branching fractions, a 1% relative uncertainty on both the signal and normalisation mode efficiencies, and a √ N uncertainty on the normalisation yield. The relative precision on R c is found to closely follow the relative precision on N (B + c → τ + ν τ ). All values shown are summarised in App. B Tab. 4. It is also possible to determine an absolute branching fraction for the signal decay, where the measured ratio of branching fractions R c is multiplied by the SM prediction for B(B + c → J/ψµ + ν µ ). In doing so, one must assume a value of |V cb | in the calculation of the normalisation branching fraction; in such a setup, an interpretation of B(B + c → τ + ν τ ) in terms of |V cb | would not be plausible. Naturally, the precision on the absolute branching fraction B(B + c → τ + ν τ ) will be impacted by the precision on the normalisation mode branching fraction prediction, which at present is 8% relative. The resulting limitation in B(B + c → τ + ν τ ) prediction can be seen in Fig. 8  (b), where the improvement with N Z is more modest than for R c due to the additional uncertainty from B(B + c → J/ψµ + ν µ ); all values shown are summarised in App. B Tab. 5. The precision of the normalisation mode branching fraction calculation is limited at present by knowledge of the B + c → J/ψ form factors, which can be improved in future measurements of the B + c → J/ψµ + ν µ decay both at LHCb and FCC-ee. In particular, complete form factor information could be determined from an angular analysis, following the approach described in Ref. [60] to measure a full set of angular coefficients. Assuming SM amplitudes in the decay would enable form factor information to derived directly from the measured angular coefficients.

Additional factors to consider
In the analysis presented above, a parametric description of the IDEA detector is employed to model key elements of the expected detector response such as momentum and impact parameter resolution. In future, studies using full simulation will be required to evaluate the complete resolution on key quantities such as the hemisphere energies, as well as to determine the reconstruction and selection efficiencies expected in genuine FCC-ee operation. Such studies should be performed under different detector design scenarios, in order to determine how aspects such as tracking, calorimetry, and vertex reconstruction may influence the expected precision on N (B + c → τ + ν τ ). The impact of particle identification (PID) should also be studied, since the current analysis is performed using combinations of three genuine charged pions only i.e. perfect PID is assumed. The pions produced in B + c → τ + ν τ decays have momenta in the 1 − 10 GeV/c range, where dE/dx, time of flight, and Cherenkov techniques can provide high discrimination between pions, kaons, and muons. Muon rejection is necessary to suppress high branching fraction semileptonic decays of beauty and charm hadrons, while kaon rejection is important since kaons are often produced in the decays of charm hadrons to multi-track final states.
It may be possible to extend the measurement to include the τ + → π + π + π − π 0ν τ mode, which also provides a 3π vertex and has a branching fraction of 50% relative to the τ + → π + π + π −ν τ decay. Due to differences in 3π kinematics and resonant structure, the selection employed in this analysis is highly inefficient on the τ + → π + π + π − π 0ν τ decay, rendering any contribution in the final sample negligible. As such, a dedicated selection using B + c → τ + ν τ with τ + → π + π + π − π 0ν τ signal MC for the MVA training and cut optimisation would be required in order to isolate these decays. Such a strategy could either include reconstruction of the additional neutral pion, or proceed with only 3π vertex reconstruction. In either case, the contribution from 3π signal decays passing the 3ππ 0 selection must be modelled.
In the analysis presented within, only combinations of pions originating from common decay vertices are considered. In real data, combinatorial background will also contribute, where one of more of the pions in the reconstructed 3π system will not originate from a common vertex. Requirements on vertex χ 2 , charged track impact parameter from the primary vertex, and track momentum, should assist in minimising such contributions. The remaining level of combinatorial background in data can be estimated using combinations of particles such π + π + π + , which are non-physical and thus represent purely random combinatorics.
Finally, samples of both inclusive and exclusive background decays should be generated at the levels expected in a dataset of size N Z ∼ 10 12 . This will enable background rejections to be accurately measured up to the required 10 10 level, with sufficient statistics remaining to model the background contributions in the signal yield fit.

Implications for New Physics
The phenomenological impact of a measurement of B(B + c → τ + ν τ ) with the precision depicted in Fig. 8 is now explored for a few NP scenarios. A measurement of the ratio with a precision of ≈ 4% is considered, and it is assumed that B(B + c → J/ψµ + ν τ ) is not affected by NP contributions, which is a well justified assumption for the models discussed below. For the external input Γ(B + c → J/ψµ + ν µ )/|V cb | 2 , we consider two benchmark scenarios: (i) a relative uncertainty of ≈ 7%, as currently obtained with LQCD form factors, and (ii) an uncertainty of ≈ 2% which could be obtained in the future by combining LQCD and experimental inputs, as discussed in Sec. 3. These two scenarios amount to an uncertainty on Γ(B + c → τ + ν τ )/|V cb | 2 of ≈ 10% and ≈ 4%, respectively, which will be considered in what follows to constrain NP contributions. 2

Effective Hamiltonian
Firstly, the most general dimension-six effective Hamiltonian encoding SM and NP contributions to the b → cτ ν τ transition is considered [37], , T } represent the effective coefficients evaluated at the renormalisation scale µ, which is taken to be µ = m b unless stated otherwise. The normalisation is such that the SM corresponds to g α = 0 for all coefficients. The B + c → τ + ν τ decay branching fraction can then be written in full generality in terms of Eq. (11), where g A(V ) ≡ g V R ∓ g V L and g P (S) ≡ g S R ∓ g S L are defined. From Eq. (12), it is clear that the B + c → τ + ν τ decays are particularly sensitive to pseudoscalar couplings of NP since they lift the helicity suppression of the SM amplitude. By considering an experimental precision of ≈ 10% on Γ(B + c → τ + ν τ )/|V cb | 2 and assuming a central value that coincides with the SM prediction, the following sensitivity (95% CL.) on the effective couplings is expected, where the couplings are assumed to be real. With an improved determination of the normalisation channel B + c → J/ψµ + ν µ , a precision of ≈ 4% on Γ(B + c → τ + ν τ )/|V cb | 2 can be assumed, which would amount to Note, in particular, that such a measurement would considerably improve the sensitivity on the coupling g P , which is only weakly constrained at present by the requirement that Γ(B + c → τ + ν τ ) should not saturate ≈ 30% of the total B + c meson width [14,15] which is determined experimentally [61]. The effective couplings defined in Eq. (11) can arise in several extensions of the SM. Of particular interest are extensions of the SM Higgs sector such as the Two-Higgs-doublet model (2HDM) [16] and specific models containing scalar and vector leptoquarks (LQs) [17,18], since they can induce the coefficients g P at low energies. The implications of the limits derived in Eq. (13) and (14) for these scenarios are now discussed.

2HDM
One of the minimal extensions of the SM consists in enlarging the Higgs sector with an additional Higgs doublet with the same quantum numbers [16]. Besides the SM-like Higgs boson, the spectrum of these models contain an extra CP -even Higgs, a neutral CP -odd scalar, as well as a charged Higgs boson that can contribute to charged-current transitions such as the one studied in this work. In order to avoid Flavor Changing Neutral Currents (FCNCs) from appearing at tree-level, one imposes via a discrete symmetry that fermions of a specific chirality and hypercharge should couple to a single Higgs doublet [62]. Four choices are then possible, which are known as 2HDM of type-I, II, X, and Y [16]. Among those, the type-II 2HDM is a popular choice which is embedded in the Minimal Supersymmetric extension of the Standard Model (MSSM), and which has a rich phenomenology in flavour-physics observables [63][64][65][66][67][68].
The tree-level contribution of the charged Higgs (H ± ) to the transition b → cτ ν τ can be matched to the effective Hamiltonian (11) generating the effective coefficients g S and g P . For the type-II 2HDM, it is found that [70] where m H ± is the H ± mass and tan β = v 2 /v 1 represents the ratio of the Higgs vacuum expectation values v 1,2 , which satisfy v SM = v 2 1 + v 2 2 = 246.2 GeV. The same expressions hold for example for the b → uτ ν τ transition by suitably replacing the quark masses.
The expected FCC-ee constraints on the plane tan β vs. m H ± are shown in Fig. 9 by using the limits on g P derived in Eq. (13) and (14) with the assumption of 10% and 4% experimental sensitivity on the value of Γ(B + c → τ + ν τ )/|V cb | 2 , respectively. In the same plot, the constraints derived from B(B → X s γ) [64] are superimposed, which set a lower limit m H ± 570 GeV at 95% CL. Furthermore, the constraints arising from the current measurement of B(B + → τ + ν τ ) [61] are shown, as well as the future prospects at Belle-II where a precision of 5% [69] is anticipated. From this plot, it is clear that the measurement of B + c → τ + ν τ at FCC-ee can probe an important region of parameter space in type-II 2HDM which will not be covered by other flavour constraints.
average of the lepton flavor universality ratios, with respect to their SM predictions, as discussed in [23] and references therein. Similar discrepancies have also been observed in the ratio R J/ψ = B(B + c → J/ψτ + ν τ )/B(B + c → J/ψµ + ν µ ) [5], but with lower statistical significance and still limited experimental precision. These discrepancies can be simultaneously explained by NP contributions to the effective coefficients defined in Eq. (11). The simplest of these explanations requires an operator with the same chirality as the SM operator, with an effective coupling within the following 1σ range, as detailed in [71] and references therein. Such an effective scenario could be induced by the tree-level exchange of the vector leptoquark U 1 = (3, 1, 2/3) or the scalar S 1 = (3, 1, 1/3), with couplings to left-handed fermions [71]. These particles are written in terms of their SM quantum numbers,  [71]. The magenta shaded region denotes the current average of R exp D /R SM D and R exp D * /R SM D * at 1σ accuracy [23]. The grey solid (dashed) lines correspond to the estimated sensitivity of 4% (10%) precision on Γ(B + c → τ + ντ )/|V cb | 2 at FCC-ee.
would imply a deviation from the SM prediction in Γ(B + c → τ + ν τ )/|V cb | 2 larger than O(10%), which could be fully probed at FCC-e, as depicted in Fig. 10.
Another viable possibility to explain the discrepancies in R D ( * ) are specific combinations of scalar and tensor operators, which are predicted in certain leptoquark models [72][73][74][75][76][77][78][79]. More precisely, the tree-level matching of the scalar leptoquarks with quantum numbers R 2 = (3, 2, 7/6) and S 1 = (3, 1, 1/3) to the effective Hamiltonian (11) can induce the combinations of couplings g S L (Λ) = +4g T (Λ) and g S L (Λ) = −4g T (Λ), respectively, at the matching scale Λ ≈ 1 TeV. After accounting for the renormalisation group running from ≈ 1 TeV to m b [80], these relations become g S L (m b ) ≈ +8.1g T (m b ) and g S L (m b ) ≈ −8.5g T (m b ), respectively, which are known to provide a good description of b → cτ ν τ data [71]. The latter scenario can explain the discrepancies in R D and R D * via purely real effective couplings, whereas the first requires purely imaginary values. The allowed ranges for these couplings are given for example in Ref. [71]. The impact of a measurement of B + c → τ + ν τ decays is even more dramatic for these scenarios due to the presence of a nonzero g S L coupling, which induces sizeable contributions to B(B c → τ + ν τ ) for the couplings required to explain R D ( * ) , as shown in Fig. 11.
Clearly, a determination of Γ(B + c → τ + ν τ )/|V cb | 2 with 4% precision can fully probe all of the leptoquark explanations of R D and R D * , as shown by the thick lines in Figs. 10 and 11. Since leptoquarks are the only viable explanations of these discrepancies which are still consistent with various flavor observables [71,81,82], high-p T LHC limits [71,[83][84][85][86][87], and electroweak precision constraints [88][89][90], the analysis presented within demonstrates the unique potential of a B + c → τ + ν τ  [71]. The magenta shaded regions denote the current experimental averages of RD and RD * at 1σ accuracy [23]. The grey region corresponds to the estimated sensitivity of 10% precision on Γ(B + c → τ + ντ )/|V cb | 2 at FCC-ee. measurement at the FCC-ee Z-pole to either confirm or refute these anomalies.

Conclusion
In this work, the prospects for a precise measurement of the branching fraction of B + c →τ + ν τ decays at the FCC-ee Z-pole are evaluated. Large simulated samples of B + c →τ + ν τ signal and associated backgrounds are generated using the FCC-ee software to simulate detector effects, reconstruct signal candidates, and perform selection optimisation and fit studies. A two-stage BDT selection is employed to reduce all sources of hadronic Z background, first using topological event-level information, and then the vertex properties of the detached τ + → π + π + π −ν τ decay to reduce the rate of b-hadron backgrounds. The sensitivities for both the branching fraction of B + c →τ + ν τ and the ratio R c = B(B + c → τ + ν τ )/B(B + c → J/ψµ + ν µ ) are estimated as a function of the number of collected Z decays, where a relative precision of around 4% is achieved for R c with N Z = 5 × 10 12 . The precision on the absolute branching fraction is limited to around 8% due to knowledge of the B + c → J/ψµ + ν µ decay form factors, which can be improved through dedicated measurements of this mode in future.
The impact of a measurement of B + c → τ + ν τ on NP scenarios is also discussed. In particular, it is shown that such a measurement at FCC-ee can constrain a large region of the (tan β, m H ± ) plane in the type-II 2HDM, which cannot be covered by other flavour-physics measurements. Recently, leptoquark models have received significant attention as the only viable explanation of the B-physics anomalies in both charged and neutral current processes. A precise measurement of the branching fraction of B + c →τ + ν τ at FCC-ee could fully probe the interpretations of R D and R D * that are permitted under existing constraints.
In summary, this work demonstrates why FCC-ee is the most well-suited environment for a measurement of the branching fraction of the B + c →τ + ν τ decay, and represents the first FCC-ee analysis to use common software tools from EDM4hep through to final analysis.  Table 2: Summary of the exclusive B-hadron background samples used for determination of the optimal BDT cuts and to model background in the signal yield fit. The yields expected for each decay mode with NZ = 5 × 10 12 are shown in the second column, and the generated sample statistics are shown in the third column. The large ratio between expected and generated statistics (fourth column) illustrates why cut-and-count efficiencies cannot be used to determine the expected background rejection achieved by the double-BDT selection; for many background modes, none of the generated events survive the optimal BDT cuts applied. The efficiency values determined using the spline parameterisation approach are given in the final column. Using splines derived from the total exclusive sample enables efficiencies down to the 10 −10 level to be evaluated.