Search for the rare decay of charmed baryon $\Lambda_c^+$ into $p \mu^+ \mu^-$ final state

A search for the nonresonant $\Lambda_c^+ \to p \mu^+ \mu^-$ decay is performed using proton-proton collision data recorded at a centre-of-mass energy of 13 TeV by the LHCb experiment, corresponding to an integrated luminosity of 5.4 fb$^{-1}$. No evidence for the decay is found in the dimuon invariant-mass regions where the expected contributions of resonances is subdominant. The upper limit on the branching fraction of the $\Lambda_c^+ \to p \mu^+ \mu^-$ decay is determined to be $2.9~(3.2) \times 10^{-8}$ at 90% (95%) confidence level. The branching fractions in the dimuon invariant-mass regions dominated by the $\eta$, $\rho$ and $\omega$ resonances are also determined.


Introduction
The flavour-changing neutral-current decay 1 Λ + c → pµ + µ − is heavily suppressed in the Standard Model (SM) as it is forbidden at tree level and is affected by the Glashow-Iliopoulos-Maiani (GIM) mechanism [1].The phenomenology of the decay is discussed in Ref. [2] and the SM predictions for q 2 differential branching fractions, where q 2 is the dimuon invariant-mass squared, are calculated in Ref. [3].The branching fraction (B) associated with the short-distance contributions to the c → uµ + µ − transition is expected to be below 10 −8 in the SM, with an enhancement up to 10 −6 due to longdistance processes involving intermediate vector-meson resonances that decay to a lepton pair [4,5].The estimation of resonant SM contributions relies on the measured values for the corresponding branching fractions and the isospin relation between ρ and ω states, to overcome the lack of experimental information on B(Λ + c → pρ).Strong phases between the resonances, and thus the effect of interference, are a priori unknown.
Interesting deviations from the SM predictions have been observed in the equivalent type of transition for the b quark, b → sℓℓ, in particular in angular observables [6].These anomalies motivate the extension of this physics program into the charm sector, which is far less explored both experimentally and theoretically.The first results, such as an angular analysis of charm mesons [7], are already available.A discussion of the discovery potential of effects beyond the SM in rare charm-baryon decays can be found in Ref. [8].
The LHCb collaboration previously searched for the Λ + c → pµ + µ − decay using the proton-proton (pp) data collected in 2011-2012 [9] and placed an upper limit of 7.7 × 10 −8 at 90% confidence level (CL), thereby improving the result from the BaBar experiment [10] by two orders of magnitude.The analysis of data collected in 2016-2018 presented in this paper follows the methodology of the previous analysis, and relies on improved techniques with the use of the CL s method [11] implementation based on unbinned likelihoods, the incorporation of background due to particle misidentification into the data model, a different definition of the short-distance signal and a more detailed treatment of the entire range of dimuon invariant mass where resonances are expected to dominate.Accordingly, the new results for the short distance contribution should not be combined with the results based on the 2011-2012 data.
The measurement is carried out in different regions of dimuon invariant mass m µµ and the Λ + c → pϕ(→ µ + µ − ) decay is used for normalisation since it has the same final state as the signal and is relatively abundant.

Detector and simulation
The LHCb detector [12,13] 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 [14], 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 [15] 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 [16].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 [17].The online event selection is performed by a trigger system [18], 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.
Simulated samples of Λ + c baryons are produced to model the effects of the detector acceptance and the imposed selection requirements.The pp collisions are generated using Pythia [19] with a specific LHCb configuration [20].Decays of unstable particles are described by EvtGen [21], in which final state radiation is generated using Photos [22].In the case of Λ + c → pµ + µ − decays, a model with final state particles equally distributed in the phase space is used.The interaction of the generated particles with the detector, and its response, is implemented using the Geant4 toolkit [23] as described in Ref. [24].The underlying pp interaction and Λ + c baryon are reused multiple times to speed up the simulation process, with independently generated Λ + c decay products [25].

Candidate selection
The selection criteria largely follow those described in Ref. [9].The Λ + c candidates are reconstructed by combining a pair of oppositely charged tracks identified as muons with an additional track identified as a proton.The three tracks are required to be compatible with originating from a common vertex and have significant χ 2  IP with respect to all PVs.Here, χ 2 IP is defined as the difference in the vertex-fit χ 2 of a given PV, reconstructed with and without the track under consideration.The Λ + c candidates must have a decay vertex displaced from any PV and be compatible with originating from one of the PVs.
The background formed by random combinations of tracks, referred to as combinatorial background, is suppressed using a multivariate classifier.A boosted decision tree (BDT) algorithm [26], as implemented in the TMVA toolkit [27,28], is trained using supervised learning with tenfold cross-validation [29] to achieve an unbiased classifier response.The BDT discrimination is done in two steps.A first BDT classifier employs topological and kinematic variables of the Λ + c candidates.A relatively loose requirement on the response of that BDT is applied as the last preselection step to enhance the purity of the normalisation data sample.After preselection, simulated samples of Λ + c → pϕ(→ µ + µ − ) decays are compared with data, and weights are determined to improve the agreement between data and simulation for the Λ + c transverse momentum and the number of tracks in the event.A second BDT classifier is trained using a set of variables that includes those of the first classifier and, in addition, kinematic, vertex and isolation variables related to the decay products.For the two classifiers, simulated events containing Λ + c → pµ + µ − decays are used as a proxy for the signal, while the background is modelled using sideband events with a pµ + µ − invariant mass within [2140,2251] and [2321,2430] MeV/c 2 .The optimisation of final requirements is performed in three dimensions: the response value of the second BDT classifier and two variables characterising the proton and muon particle identification (PID).In this optimisation, the simulation is corrected for imperfections in modeling of the PID performance [30].The final requirements are determined using pseudoexperiments by finding the best expected upper limit (using the CL s method) on the ratio of branching fractions of the signal relative to the normalisation channel.At this stage, only two background sources are expected to significantly contribute: combinatorial background and Λ + c → pπ + π − decays where both pions are misidentified as muons.In each pseudoexperiment, pseudodata for the signal and these two background components are drawn from the three-dimensional probability density functions determined after the preselection stage.
Six dimuon invariant-mass ranges are considered as presented in Table 1.The nonresonant signal decay is searched for in two m µµ regions, low-m and high-m, where the expected contributions of resonances are subdominant.In the following, the combination of these two ranges is referred to as the signal region.The non signal range is divided into regions where either the η, ρ, ω or ϕ resonance dominates.The branching fraction results quoted for a given resonant region are labeled by the resonant decay expected to dominate in that region.However, due to the relatively low sample size, no attempt to separate the resonant contributions within a given dimuon invariant-mass range is made.Therefore it has to be emphasised that the results must be interpreted as corresponding to that particular m µµ region rather than to a specific resonant decay.

Normalisation and invariant-mass fit
The Λ + c → pµ + µ − branching fraction is determined from the measured signal yield N sig according to where the normalisation factor α is given by Here, N norm is the yield of Λ + c → pϕ(→ µ + µ − ) decays, while ϵ sig and ϵ norm represent the corresponding total efficiencies for signal and normalisation channels, respectively.In the case of the ratio of branching fractions B(Λ , the corresponding normalisation factor α takes the form ϵ norm /(ϵ sig × N norm ).The ratio of efficiencies ϵ norm /ϵ sig for the signal regions is compatible with unity within statistical uncertainty.
An extended unbinned maximum-likelihood fit to the pµ + µ − invariant mass is performed in each m µµ region, where the signal shape of Λ + c candidates is described by a Crystal Ball function [31].The background is described using two distinct components: the combinatorial background, modeled by an exponential shape, and the Λ + c → pπ + π − peaking background, represented by the sum of a Crystal Ball function and a Gaussian function The shapes of both the signal and the Λ + c → pπ + π − background are determined from simulation for each m µµ region and fixed in the fit to data.Weights are applied to simulated Λ + c → pπ + π − events to match the phase-space distribution observed in Ref. [32].In each m µµ region, there are four parameters in the fit.Three of them are free parameters: the signal yield, the combinatorial background yield and the exponent describing the shape of the combinatorial background.The fourth one corresponds to the Λ + c → pπ + π − background.Its total yield in all regions is a free parameter, and the fraction in each region is constrained based on simulation.
The results of the fit to the data for the signal region and for the low-m and high-m regions are presented in Fig. 1.The results of the fit for the normalisation mode and other regions dominated by resonances are presented in Fig. 2. The yields and the Λ + c signal significance for the six m µµ regions are reported in Table 2, where the significance is determined using Wilks' theorem [33].For the signal regions the upper limits for the corresponding branching fractions are determined.The residual contribution from resonant decays into signal regions is estimated from simulation neglecting interference effects.For the high-m region, where the significance is the highest, only 2.8 events are attributed to Λ + c → pϕ(→ µ + µ − ) decays while for the low-m region the expected contribution is negligible.

Systematic uncertainties
The systematic uncertainties taken into account in the determination of the upper limit on the branching fraction are summarised in Table 3.The two main contributions are the statistical uncertainty of the simulation samples used to determine the efficiency ratio and the uncertainty on the yield of the normalisation mode.The latter is dominated by the statistical uncertainty on the yield returned from the fit.Other sources of systematic uncertainty on the efficiency ratio are associated with the following: • The level of agreement between data and simulation.The maximum difference between ratios of efficiencies determined for a few variants of the baseline procedure of weighting simulation to match the data is taken as a systematic uncertainty.
• The efficiency determination from simulation.The difference in the efficiency obtained including or excluding Λ + c candidates for which one or more decay products are not correctly matched to the generated particles, is taken as a systematic uncertainty.The procedure is applied consistently for the signal and normalisation modes.
• The correction of PID variables in simulation.The difference of the efficiency ratio for the baseline and an alternative parameterisation of the PID calibration procedure [30] is taken as a systematic uncertainty.
These three sources added in quadrature amount to a relative uncertainty of 1.7%.The systematic uncertainties related to the fit model are estimated using pseudoexperiments.An alternative model for each component is used to generate pseudodata, which are then fitted using the baseline model described in Sec. 4. A Johnson function [34] is used as an alternative shape for the signal, while the Λ + c → pπ + π − background is parameterised as the sum of a Johnson and a Gaussian function, and a linear function is used as an alternative model for the combinatorial background.The difference between the generated signal yield and the average fitted yield is taken as systematic uncertainty.In addition, the uncertainty related to the Λ + c → pπ + π − decay model is determined by applying different corrections to the Λ + c → pπ + π − simulated sample and determining an alternative invariant-mass shape.The bias due to maximum-likelihood fit is determined from the same data model used for pseudodata generation and fitting.The signal yield difference between the generated value and the average of fitted values from the pseudoexperiments is taken as a systematic uncertainty.
The systematic uncertainties associated with the ratio of branching fractions in the resonant regions are determined in similar way as for the signal regions and are presented in Table 3.

Signal regions
No significant contribution of Λ + c → pµ + µ − decays is observed in either of the two signal regions.Upper limits on the branching fractions and their ratios are determined using the CL s method implemented in the GammaCombo package [35].The systematic uncertainties are included in the determination of CL s .The following upper limit in the signal region is determined to be Using the values of the branching fractions for Λ + c → pϕ and ϕ → µ + µ − decays from Ref. [36] and including their uncertainties in the CL s determination, an upper limit on the branching fraction is determined to be B(Λ + c → pµ + µ − ) < 2.9 (3.2) × 10 −8 at 90% (95%) CL.
The upper limit placed using data from the high-m region is driven by the signal yield whose significance is close to but lower than 3σ, while in the low-m region the invariantmass distribution is fully compatible with no signal.Assuming a decay model with final state particles equally distributed in the phase space, the values of upper limits extrapolated from the signal region to the full phase space are determined to be 2) × 10 −8 at 90% (95%) CL.
In spite of the enlarged data sample, the extrapolated upper limit on the branching fraction is only slightly better than that determined for the 2011-2012 data [9], 7.7 (9.6) × 10 −8 at 90% (95%) CL.This relatively small improvement is driven mainly by the 2.8σ effect observed in the high-m region in this analysis.

Resonant regions
Excluding the normalisation channel, the significance of the Λ + c signal exceeds the 5σ threshold in the ρ and ω regions while the significance for the Λ + c → pη(→ µ + µ − ) decay is at 3σ.The corresponding branching fractions are determined using Eqs.1,2 where the decay parameters for the nonresonant decays are replaced by those for the resonant decays.It has to be emphasised that disentangling the resonant components to account for interference effects has not been attempted.The efficiency for a given m µµ region is determined using simulation samples of the corresponding resonant decay, thus neglecting possible contributions from other resonances.Under this assumption, the ratio of branching fractions with respect to the normalisation channel are determined for the ω, ρ and η regions to be = 0.240 ± 0.030 (stat.)± 0.018 (syst.), = 0.229 ± 0.051 (stat.)± 0.022 (syst.), = 0.032 ± 0.013 (stat.)± 0.004 (syst.).
The corresponding branching fractions are determined to be where the second uncertainties are induced by the ratio of branching fractions presented in Table 3 and the third uncertainties are due to the following branching fractions: [36], which were used as external inputs.

Summary
A search for the nonresonant Λ + c → pµ + µ − decay is presented using data samples of proton-proton collisions collected by the LHCb experiment in 2016-2018, corresponding to an integrated luminosity of 5.4 fb −1 .
In a narrow region around the ω resonance, a Λ + c signal is observed with a significance of more than 7σ, corresponding to a branching fraction (9.82 ± 1.23 (stat.)± 0.73 (syst.)± 2.79 (ext.))× 10 −4 , compatible with and more precise than the result of the previous analysis for the data collected in 2011-2012 [9].
Excluding the region around the ω meson, a signal in the ρ region is found with a significance of 5.6σ.The corresponding branching fraction determined with respect to the normalisation channel is measured to be (1.52 ± 0.34 (stat.)± 0.14 (syst.)± 0.24 (ext.))× 10 −3 , giving strong evidence for the + c → pρ decay, qualified by the fact that interference effects are not accounted for.For the η region the branching fraction was measured to be (1.67 ± 0.69 (stat.)± 0.23 (syst.)± 0.34 (ext.))× 10 −3 with a significance of 3.0σ.The value is compatible with the world average [36]

Figure 1 :
Figure 1: Distributions of pµ + µ − invariant mass for (top) signal, (bottom left) low-m and (bottom right) high-m regions.The results of the fit are overlaid.

Table 2 :
Results of the fit to the pµ + µ − invariant-mass distribution for various m µµ regions.The significance of the Λ + c → pµ + µ − component is also given.

Table 3 :
Summary of systematic uncertainties for the signal regions and for the η, ρ and ω resonant regions.