Search for $B_c^+\to\pi^+\mu^+\mu^-$ decays and measurement of the branching fraction ratio ${\cal B}(B_c^+\to\psi(2S)\pi^+)/{\cal B}(B_c^+\to J/\psi \pi^+)$

The first search for nonresonant $B_c^+\to\pi^+\mu^+\mu^-$ decays is reported. The analysis uses proton-proton collision data collected with the LHCb detector between 2011 and 2018, corresponding to an integrated luminosity of 9 fb$^{-1}$. No evidence for an excess of signal events over background is observed and an upper limit is set on the branching fraction ratio ${\cal B}(B_c^+\to\pi^+\mu^+\mu^-)/{\cal B}(B_c^+\to J/\psi \pi^+)<2.1\times 10^{-4}$ at $90\%$ confidence level. Additionally, an updated measurement of the ratio of the $B_c^+\to\psi(2S)\pi^+$ and $B_c^+\to J/\psi \pi^+$ branching fractions is reported. The ratio ${\cal B}(B_c^+\to\psi(2S)\pi^+)/{\cal B}(B_c^+\to J/\psi \pi^+)$ is measured to be $0.254\pm 0.018 \pm 0.003 \pm 0.005$, where the first uncertainty is statistical, the second systematic, and the third is due to the uncertainties on the branching fractions of the leptonic $J/\psi$ and $\psi(2S)$ decays. This measurement is the most precise to date and is consistent with previous LHCb results.


Introduction
The B + c meson is made of the heaviest quarks that bind to form hadrons in the Standard Model (SM): beauty and charm quarks.The presence of two heavy quarks distinguishes the system from other B mesons in theoretical calculations.Hence, measurements of the Table 1: Ranges of dimuon invariant mass or mass squared (q 2 ) used in the analysis.
For each of the intervals, a fit is performed using the B + c candidate invariant mass, m(π + µ + µ − ), as discriminating observable.The yield, relative to that for the B + c → J/ψπ + normalisation mode, is converted to a result on the relative branching fraction of the given q 2 bin, where N indicates a yield, ε indicates the efficiency determined from simulation with data-driven corrections, and B(J/ψ → µ + µ − ) is the known branching fraction of the J/ψ → µ + µ − decay [21].Here, B(B + c → π + µ + µ − ) indicates the B + c → π + µ + µ − differential branching fraction integrated over the q 2 range relevant for the given bin.

Detector and simulation
The LHCb detector [22,23] 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 pp interaction region [24], a large-area silicon-strip detector located upstream of a dipole magnet with a bending power of about 4 T m, and three stations of silicon-strip detectors and straw drift tubes [25,26] 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.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.Different types of charged hadrons are distinguished using information from two ring-imaging Cherenkov detectors [27].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 [28].The online event selection is performed by a trigger [29,30], which consists of a hardware stage, based on information from the calorimeter and muon systems, followed by a two-level software stage, which reconstructs the full event.
Simulation is used to optimise the event selection procedure, to model the shape of the B + c candidate invariant-mass distributions and to estimate efficiencies accounting for the effects of the detector acceptance, reconstruction and selection criteria.In the simulation, pp collisions are generated using Pythia [31] with a specific LHCb configuration [32].The production of B + c mesons is simulated using the dedicated generator BcVegPy [33].Decays of unstable particles are described by EvtGen [34], in which final-state radiation is generated using Photos [35].The interaction of the generated particles with the detector, and its response, are implemented using the Geant4 toolkit [36] as described in Ref. [37].The underlying pp interaction is reused multiple times, with an independently generated signal decay for each [38].
The B + c candidates reconstructed in simulation are weighted to correct for discrepancies between data and simulation associated with the particle-identification [39], trackreconstruction [40] and hardware trigger [41] efficiencies.The simulation is also corrected such that the B + c lifetime corresponds to the current experimental value [21,42,43].Additional corrections are applied to account for discrepancies in B + c production kinematics, event track multiplicity and other observables used in the selection of B + c candidates.These corrections are obtained using a multivariate algorithm [44], which is trained using B + c → J/ψπ + decays in background-subtracted data and simulation.After the corrections are applied, the simulated distributions of all variables used in the analysis are in good agreement with the data.

Candidate selection and background sources
Candidate B + c → π + µ + µ − decays are triggered in an identical manner as described in Ref. [45] for B + decays to the same final states.The hardware stage of the trigger selects events containing at least one muon with high p T .The software stage of the trigger selects events containing at least one high-p T muon that is inconsistent with originating from any PV.The events must contain at least one secondary vertex (formed by two or more of the final-state particles) that is significantly displaced from every PV.A multivariate algorithm [46,47] is used to identify secondary vertices consistent with the decay of a b hadron.
The initial stages of the offline selection are also similar to those for B + → π + µ + µ − decays [45], except that less stringent impact parameter and flight distance requirements are imposed to account for the shorter B + c lifetime compared to that of the B + meson.The B + c candidates are formed from pairs of well-reconstructed oppositely charged tracks that are identified as muons together with a track identified as a pion, and are required to have a good-quality vertex.Each B + c candidate is associated with the PV giving the smallest value of χ 2 IP , which is defined as the difference in the vertex-fit χ 2 of the PV reconstructed with and without the candidate.Each B + c candidate must be consistent −3 39 ± 7 11.0 < q 2 < 12.5 GeV 2 −18 +7 −10 model includes two components: signal decays and combinatorial background.The signal model is validated using simulation and corresponds to the sum of two Gaussian functions, one with power-law tails on both sides of the distribution [52].The tail parameters are fixed to the values obtained from simulation.The signal model includes a global shift of peak position and a global scaling factor for the width of the distribution, relative to the values found in simulation.The peak offset and width scale factor are obtained from a fit to the resonant B + c → J/ψπ + decays in data.The combinatorial background is described by an exponential function, whose exponent is allowed to vary in the fit to data.
Figure 1 shows the B + c candidate invariant-mass distributions of selected B + c → π + µ + µ − candidates, with results of the fits superimposed.For the intermediate q 2 interval, the fit favours a negative signal yield due to a lack of candidates in the region close to the B + c mass peak position and a feature of extended maximum-likelihood fits that prefers negative yields in such situations [53].Table 2 summarises the yields obtained from the fits.
For the resonant B + c → J/ψπ + decay, the fit model includes four components: signal B + c → J/ψπ + decays, misidentified B + c → J/ψK + decays, partially reconstructed background from B + c → J/ψρ + decays and combinatorial background.The signal, misidentified and partially reconstructed backgrounds are each parameterised by the sum of two Gaussian functions, one with power-law tails.The tail parameters of each distribution are fixed from simulation.The peak position and width of the distributions are allowed to vary in the fit to the data by a global offset and scale factor that is shared between the components.The combinatorial background model is an exponential function, whose exponent is allowed to vary in the fit to data.In total, the B + c → J/ψπ + fit includes seven free parameters: the yields of the four components, the global peak position shift and width scaling factor, and the exponent of the combinatorial background.
Figure 2 shows the distributions of selected B + c → J/ψπ + candidates for both WPs, with fit projections overlaid.Table 3 summarises the yields obtained from the fits.The ratio of B + c → J/ψπ + yields at the two WPs is found to be consistent with the expectation based on the efficiencies in simulation.The yields of the misidentified B + c → J/ψK + component at the two WPs are found to be consistent with the expectations based on the misidentification rates in simulation and the measured ratio of branching fractions B(B + c → J/ψK + )/B(B + c → J/ψπ + ) [51].Similarly, for the resonant B + c → ψ(2S)π + decay, the fit includes four components:

fb
All Figure 1: Reconstructed π + µ + µ − invariant-mass distributions for the selected B + c → π + µ + µ − candidates for each q 2 interval and for all intervals combined, with results of the fit described in the text overlaid.
background and combinatorial background.The analytical functions of the fit models used for the B + c → ψ(2S)π + fit are the same as those used for the B + c → J/ψπ + fit, but with different parameters.The signal model, and the models for the misidentified and the partially reconstructed background have tail parameters which are fixed to values obtained from simulation, and they also include a global shift of peak position and a global scaling factor for the widths of the distributions.However, the global peak position shift and width scaling factor are constrained in the fit to data to be consistent with the values obtained from the B + c → J/ψπ + fit.Thus, the

fb
Figure 2: Reconstructed π + µ + µ − invariant-mass distributions for the selected B + c → J/ψπ + candidates at the (left) π + µ + µ − and (right) ψ(2S)π + working points, with the results of the fits overlaid.For the right figure, the π + µ + µ − invariant mass is calculated after constraining the dimuon mass to the known J/ψ mass, thereby improving the signal resolution.B + c → ψ(2S)π + fit includes five free parameters: the yields for the four components and the exponent of the combinatorial background.The B + c → ψ(2S)π + fit is performed to the B + c mass distribution after constraining the dimuon mass to the known ψ(2S) mass.Figure 3 shows the distribution of selected B + c → ψ(2S)π + candidates, with fit projections overlaid.Table 4 summarises the yields obtained from the fit.

Component
Yield Table 5: Efficiency ratios between signal and normalisation modes.The values of ε J/ψπ + /ε π + µ + µ − are provided for each q 2 interval and for all intervals combined.The uncertainty on the ratio combines statistical and systematic effects.
Signal decay mode

Efficiencies and systematic uncertainties
The efficiency ratios between signal and normalisation modes are obtained from simulation accounting for the geometrical acceptance of the detector as well as effects related to the triggering, reconstruction and selection of the B + c candidates.Table 5 lists the efficiency ratios between signal and normalisation modes.For the nonresonant signal mode, the efficiency denominator includes only the generated events in the respective q 2 interval.The uncertainties on the efficiency ratios take into account the simulation sample size, uncertainties on the weights applied to the simulation, the matching between reconstructed and generated particles in the simulation, variations of the software trigger requirements, and the uncertainty on the known B + c lifetime.All variations are made consistently for the signal and normalisation modes to avoid overestimation of the uncertainty on the efficiency ratio.For the nonresonant signal mode, the impact of the signal decay model assumed in the simulation is additionally considered.
The systematic uncertainties associated with the weights are evaluated by varying all weights within their uncertainties and by varying the binning scheme used to estimate them.The systematic uncertainty associated with the multivariate weighting algorithm (see Sec. 2) is evaluated by comparing the results obtained with the default and with an alternative algorithm.The default algorithm is trained to correct for discrepancies between data and simulation associated with the event track multiplicity and with the p T and the vertex quality of the B + c candidates.The alternative algorithm is trained using the IPs of the two muons as additional inputs.
The systematic uncertainty associated with the matching between reconstructed and generated particles in the simulation is evaluated by comparing the efficiencies obtained including or excluding B candidates for which one or more decay products are not correctly matched.The systematic uncertainty associated with variations of the software trigger requirements that are not mirrored by the simulation is evaluated by comparing the efficiencies obtained by applying the tightest thresholds and by applying average thresholds within each data-taking period.The systematic uncertainty associated with the B + c lifetime is evaluated by varying the B + c lifetime in simulation within its uncertainties [21].The nonresonant signal decays are simulated with a phase-space distribution of the finalstate particles.However, the true q 2 distribution of B + c → π + µ + µ − decays is unknown, and there is no clear theoretical guidance for its expectation.A systematic uncertainty is assigned due to the change of efficiency over q 2 , shown in Fig. 4. The uncertainty is taken to be the root-mean-square variation of the distribution of the efficiency values in narrow q 2 bins within each interval.The approach of determining the differential branching fraction within q 2 bins helps to reduce this uncertainty.
A similar uncertainty arises due to the different possible polarisations in the dimuon system, and is illustrated by the different coloured distributions in Fig. 4. The two extreme possibilities correspond to the dimuon system forming a scalar state, which is unpolarised, and a vector state that is longitudinally polarised.The two cases are characterised by different cos θ l distributions, where the helicity angle θ l is defined as the angle between the µ + momentum and the opposite of the B + momentum in the dimuon frame: the scalar dimuon state corresponds to a flat cos θ l distribution, while the vector dimuon state corresponds to a 3 4 sin 2 θ l distribution.The difference in efficiency between the two extreme cases is taken as the associated systematic uncertainty.
For the resonant modes, the effect of the multivariate weighting algorithm dominates the systematic uncertainty on the efficiency ratio, while for the nonresonant signal mode, the systematic uncertainty associated with efficiency variation within the q 2 intervals dominates the uncertainty on the efficiency ratio.The remaining systematic uncertainties cancel out almost fully in the determination of the efficiency ratios.For all measured quantities, the systematic uncertainties are smaller than the statistical uncertainties.
The yields obtained from the fits described in the previous section can be affected by the fit model choice and by the assumption of the polarisation of the partially reconstructed backgrounds.To study the effect of the fit model choice, each fit is performed in two configurations: using the baseline fit model and using an alternative fit model.In the alternative fit model, the analytical function used for the signal and the misidentified background models is replaced with the sum of two Gaussian functions, one with a power-law tail on the left side of the distribution, and the other with a power-law tail on the right.The analytical function for the partially reconstructed backgrounds is also replaced, but with the sum of two Gaussian functions.In all cases the difference between the results obtained with the baseline and alternative models is found to be negligible.
For the fits to nonresonant signal candidates, the possible contribution from partially reconstructed B + c → ρ + µ + µ − decays is found to be negligible by performing fits to sideband data of background-only models including a B + c → ρ + µ + µ − background component.The sideband data includes B + c candidates in the fit range excluding those in the signal region (see Sec. 3).For the resonant modes, the ρ + meson is assumed to be unpolarised.However, the polarisation of the ρ + meson can affect the momentum of the missing pion and hence the B + c candidate mass shape of the partially reconstructed backgrounds.To study the effect of the ρ + polarisation, the fits are repeated assuming either full longitudinal or full transverse ρ + polarisation.The difference in the results for the two Figure 4: Efficiency ratio between the B + c → J/ψπ + and nonresonant B + c → π + µ + µ − decay modes as a function of q 2 , evaluated for the two extreme possibilities of the dimuon system forming a scalar state, which is unpolarised, and a vector state that is longitudinally polarised.The shaded q 2 intervals, which contain the contributions from the J/ψ and ψ(2S) resonances, are not used in the analysis of nonresonant decays.
configurations is found to be negligible.
Upper limits on the branching fraction ratios are obtained following the Feldmann-Cousins prescription [54]: pseudoexperiments are generated for various values of R π + µ + µ − /J/ψπ + and the resulting distribution of measured R π + µ + µ − /J/ψπ + is used to form confidence belts.Figure 5 shows confidence belts at 90% and 95% confidence level (CL).Table 6 gives the results for R π + µ + µ − /J/ψπ + and the obtained limits.To assess the impact of the systematic uncertainties, the fits are repeated fixing the nuisance parameters to their central values.Figure 6 summarises the obtained limits on the normalised differential branching fraction.As further checks the procedure is repeated restricting the signal yield to positive values, or performing nonextended maximum-likelihood fits.No significant changes in the obtained upper limits are found.

Summary
A search for nonresonant B + c → π + µ + µ − decays is performed together with an updated measurement of the ratio of the B + c → ψ(2S)π + and B + c → J/ψπ + branching fractions.The analysis uses proton-proton collision data collected with the LHCb detector between 2011 and 2018, corresponding to an integrated luminosity of 9 fb −1 .No evidence for an excess of signal events over background is observed for nonresonant B + c → π + µ + µ − decays and an upper limit is set on the branching fraction ratio at 90% confidence level.This is the first limit on B + c decays mediated only by annihilation diagrams into a semileptonic final state.For the resonant B + c → ψ(2S)π + mode, the branching fraction ratio is measured to be

All
Figure 5: Confidence belts generated using pseudoexperiments according to the Feldman-Cousins prescription [54] for each q 2 interval and for all intervals combined.The vertical dashed line shows the central value from the fit to data.
where the third uncertainty is due to limited precision of the known branching fractions for J/ψ and ψ(2S) leptonic decays [21].This measurement is consistent with, and supersedes, previous LHCb results on the same quantity [19,20] and is the most precise to date. Figure 6: Upper limits on the normalised differential branching fraction for nonresonant B + c → π + µ + µ − decays as a function of q 2 .The solid lines show the results for each q 2 bin, while the dashed lines show the results for all bins combined.

Table 2 :
Yields for the signal nonresonant B + c → π + µ + µ − decay (N π + µ + µ − ) and combinatorial background (N comb ) obtained from the fits to data described in the text, with statistical uncertainties only.

Table 3 :
Yields obtained from the B + c → J/ψπ + fits to data at the two working points as described in the text, with statistical uncertainties only.

Table 4 :
Yields obtained from the B + c → ψ(2S)π + fit to data described in the text, with statistical uncertainties only.

Table 6 :
Results for R π + µ + µ − /J/ψπ + , where the first uncertainties are statistical and the second are systematic, together with upper limits (ULs) at 90% and 95% CL.