Study of charmonium and charmonium-like contributions in $B^+ \rightarrow J/\psi \eta K^+$ decays

A study of $B^+ \rightarrow J/\psi \eta K^+$ decays, followed by $J/\psi \rightarrow \mu^+ \mu^-$ and $\eta \rightarrow \gamma \gamma$, is performed using a dataset collected with the LHCb detector in proton-proton collisions at centre-of-mass energies of 7, 8 and 13 TeV, corresponding to an integrated luminosity of 9 fb$^{-1}$. The $J/\psi\eta$ mass spectrum is investigated for contributions from charmonia and charmonium-like states. Evidence is found for the $B^+\rightarrow \left( \psi_2(3823) \rightarrow J/\psi \eta \right) K^+$ and $B^+\rightarrow \left( \psi(4040) \rightarrow J/\psi \eta \right) K^+$ decays with significance of 3.4 and 4.7~standard deviations, respectively. This constitutes the~first~evidence for the $\psi_2(3823) \rightarrow J/\psi \eta$ decay.

The J/ψη final state is well suited to search for the hypothetical X C state. Searches for the X C → J/ψη decay have been performed by the BaBar and Belle collaborations using B + → J/ψηK + decays [66,67]. No X C → J/ψη signal is observed, and upper limits at 90% confidence level (CL) on the product of branching fractions for the B + → X C K + and X C → J/ψη decays are set to be 7.7 × 10 −6 (BaBar) and 3.8 × 10 −6 (Belle).

Event selection
Candidate B + → J/ψηK + decays are reconstructed through the J/ψ → µ + µ − and η → γγ decay modes. A loose initial selection is applied to reduce the background. The criteria are chosen to be similar to those used in previous LHCb studies [21, 94-98]. Subsequently, a multivariate estimator based on an artificial neural network algorithm [99,100], configured with a cross-entropy cost estimator [101], in the following referred to as the MLP classifier, is applied. Muon and kaon candidates are identified by combining information from the Cherenkov detectors, calorimeters and muon detectors [102] associated to the reconstructed tracks. Transverse momenta of muon candidates are required to be greater than 550 MeV/c. To reduce combinatorial background only tracks that are inconsistent with originating from any reconstructed PV in the event are considered. Pairs of oppositely charged muons consistent with originating from a common vertex are combined to form J/ψ → µ + µ − candidates. The reconstructed mass of the pair is required to be between 3.056 and 3.136 GeV/c 2 .
Photons are reconstructed from clusters in the electromagnetic calorimeter that have transverse energy larger than 500 MeV and are not associated with reconstructed tracks [85,103]. Photon identification is based on the combined information from electromagnetic and hadronic calorimeters, scintillation pad and preshower detectors and the tracking system. Candidate η → γγ decays are reconstructed as diphoton combinations with mass within ±60 MeV/c 2 of the known η mass [34] and transverse momentum greater than 1.5 GeV/c.
The selected J/ψ candidates are combined with K + and η candidates to form B + candidates. Each B + candidate is associated with the PV that yields the smallest χ 2 IP , where χ 2 IP is defined as the difference in the vertex-fit χ 2 of a given PV reconstructed with and without the charged tracks that form the B + candidate under consideration. To improve the B + mass resolution a kinematic fit [104] is performed. This fit constrains the µ + µ − and γγ masses to the known J/ψ and η mesons masses [34], respectively, and the B + candidate to originate from its associated PV. The proper decay time of the B + candidate is required to be greater than 200 µm/c to suppress the large combinatorial background.
A further selection based on the MLP classifier reduces the combinatorial background to a low level whilst retaining a high signal efficiency. Variables included in the classifier are related to the reconstruction quality, kinematics and decay time of the B + candidates, kinematics of the final-state particles and a variable that characterises kaon identification. The classifier is trained using simulated samples of B + → J/ψηK + decays as signal proxy. The B + → J/ψηK + candidates from data with mass, m J/ψηK + , ranging between 5.4 and 5.7 GeV/c 2 , are used as background proxy. To avoid introducing a bias in the MLP evaluation, a k-fold cross-validation technique [105] with k = 13 is used.
The requirement on the MLP classifier is chosen to maximize the figure-of-merit defined as S/ √ B + S, where S represents the expected signal yield from simulation, and B stands for the background yield obtained by fitting the data. The expected signal yield is estimated as S = εS 0 , where ε is the efficiency of the requirement on the MLP classifier determined from simulation, and S 0 is the signal yield obtained from the fit to the data when no requirement is applied. The mass distribution of selected B + → J/ψηK + candidates is shown in Fig. 1.

B + → J/ψηK + signal and J/ψη mass spectrum
The B + → J/ψηK + signal yield is determined using an extended unbinned maximum-likelihood fit to the J/ψηK + mass distribution with signal and background components. The signal is modelled by a modified Gaussian function with power-law tails on both sides of the distribution [106,107], referred to hereafter as F S . The tail parameters of the modified Gaussian function are fixed from simulation, while the peak position and resolution are allowed to vary in the fit. The combinatorial background is parametrised with an exponential function. The fit result is overlaid in Fig. 1 and the signal yield is found to be N B + →J/ψηK + = (5.39 ± 0.16) × 10 3 .
The search for the B + → (X → J/ψη) K + signal is performed using extended unbinned maximum-likelihood fits to the background-subtracted J/ψη mass spectrum. The sP lot technique [108] is used for the background subtraction using the J/ψηK + mass as the discriminating variable. To improve the J/ψη mass resolution and to eliminate a small correlation between the m J/ψηK + and m J/ψη variables, the J/ψη mass is computed using a kinematic fit [104] that constrains the mass of the B + candidate to its known value [34]. For easier parametrisation of the nonresonant component, the fit to the J/ψη mass distribution is performed separately in four different overlapping mass regions. For X masses below 3.875 GeV/c 2 , a fit of the lowestmass region, 3.65 < m J/ψη < 3.90 GeV/c 2 , is performed.
The fit function to the lowest-mass region consists of three components: 1. the decay of interest B + → (X → J/ψη) K + , referred to as C X component; 3. the B + → (J/ψη) NR K + decays with no resonances in the J/ψη system, and referred to as C NR component.
The C ψ(2S) component and the C X contribution for X states with the natural width negligible with respect to the detector resolution (referred as narrow) are modelled using the F S shape. The tail parameters of all the F S functions are fixed from simulation, while the peak position and resolution parameter for the C ψ(2S) component are allowed to vary in the fit. The ratio of the resolution parameters for the C X and C ψ(2S) components is fixed at the value obtained from simulation. This procedure also accounts for a small imperfection in the modelling of the J/ψη mass resolution in the simulation [29, 30, 33].
The nonresonant component C NR is parameterised with a product of the phase-space function describing a two-body system out of the three-body final state [109] and a positive first-order polynomial function. For X masses above 3.875 GeV/c 2 , the fit is performed simultaneously in two J/ψη mass regions, one containing the X mass and the other the ψ(2S) state. For the X mass region the fit function consists of a pair of components, C X and C NR , while the ψ(2S) region fit includes the C ψ(2S) and C NR components. The parameters for two C NR components are independent in the two regions.
No significant signal is found for B + → (X → J/ψη) K + decays occurring via a hypothetical narrow X particle in the 3.7 < m X < 4.7 GeV/c 2 mass region. Fit results without the C X component are shown in Fig. 2, illustrating that no sizeable contribution from decays with a narrow intermediate X state is required to describe the data. To quantify the absence of the B + → (X → J/ψη) K + signal, fits are performed with the mass m X of the hypothetical X particle fixed to a value that is scanned across the whole available J/ψη mass range. The yield of the C X component, N X , is parametrised using the yield of the C ψ(2S) component, N ψ(2S) , and the ratio of branching fractions, F X , defined by Eq. (1), as where R ε is the ratio of total efficiencies for the B + → (X → J/ψη) K + and B + → (ψ(2S) → J/ψη) K + channels, described in Sec. 5. The parameters N ψ(2S) and F X are left to vary in the fit, and the uncertainty on the ratio R is included in the fit through a Gaussian constraint. A second set of fits exploits an alternative parametrisation that allows for the determination of the product of the branching fractions B X , defined by Eq.
In addition to the search for decays with a narrow hypothetical X state, a search is performed for decays mediated by known conventional charmonium or charmonium-like states, including the hypothetical X C state and the neutral partner of the Z c (4430) + state. For the latter it is assumed that the mass and width are the same as for its charged partner [78][79][80][81][82], while for the X C state the mass and width are assumed to be the same as for the χ c1 (3872) state [29,30]. The C X component is parametrised with a relativistic S-wave Breit-Wigner shape convolved with the F S function. For each resonance with a mass larger than 3.9 GeV/c 2 a fit range is chosen individually depending on the resonance mass and width. The uncertainties on the resonance parameters are included in the fits using Gaussian constraints. For the ψ 2 (3823) state, where only the upper limit for the natural width is known [30], a natural width of 1 MeV is assumed. The background-subtracted J/ψη mass spectra in the corresponding ranges, together with the fit results, are shown in Figs. 3 to 5. For the B + → (ψ 2 (3823) → J/ψη) K + and B + → (ψ(4040) → J/ψη) K + decays, signals with a statistical significance of 3.4 and 9.0 standard deviations, respectively, are seen. No evidence for other decays is obtained.

Efficiency and systematic uncertainty
For each considered value of m X , the efficiency ratio R ε from Eq. (3a) is calculated as  where the total efficiency ε for each decay is calculated from the product of the detector acceptance, the reconstruction and selection efficiencies for decays within the detector acceptance, and the trigger efficiency for decays passing the selection criteria. All efficiencies are calculated using simulation, as described in Sec. 2. The finite size of the simulation samples contributes to the uncertainty on the R ε (m X ) ratio. Since signal and normalisation decays share the same final state, many systematic uncertainties cancel in the ratio R ε . The remaining nonnegligible uncertainties are listed in Table 1. A large class of systematic uncertainties is associated to the corrections applied to the simulation. The finite size of the B + → J/ψK + signal sample used for correction of the simulated, p T and y spectra of B + mesons, induces a corresponding uncertainty. In turn, the variation within their uncertainties induces small changes in the ratio R ε . The corresponding spread of these changes amounts to 0.1% and is taken as systematic Table 1: Relative systematic uncertainties for the efficiency ratio R ε . When an uncertainty is found to be dependent on the J/ψη mass, the corresponding range is shown. The total uncertainty is obtained as the quadratic sum of the individual contributions.

Source
Uncertainty The kaon identification variable used for the MLP estimator is drawn from calibration data samples accounting for the dependence on particle kinematics and track multiplicity. Systematic uncertainties in this procedure arise from the limited statistics of both the simulation and calibration samples, and the modelling of the identification variable. The limitations due to both simulation and calibration sample size are evaluated by bootstrapping to create multiple samples, and repeating the procedure for each sample. The impact of potential mismodelling of the kaon identification variable is evaluated by describing the corresponding distributions using density estimates with different kernel widths [92]. For each of these cases, alternative efficiency maps are produced to determine the associated uncertainties. A systematic uncertainty of 1% is assigned from the observed differences with alternative efficiency maps.
There are residual differences in the reconstruction efficiency of charged-particle tracks that do not cancel completely in the ratio R ε due to the slightly different kinematic distributions of the final-state particles. The track-finding efficiencies obtained from simulated samples are corrected using calibration channels [93]. The uncertainties related to the efficiency correction factors are propagated to the ratios of the total efficiencies using pseudoexperiments and are found to be less than 0.15% for the considered values of the m X parameter. Differences between data and simulation of photon reconstruction efficiencies are studied using a large sample of B + → J/ψ (K * + → K + (π 0 → γγ)) decays [94,95,110]. The associated systematic uncertainty largely cancels.
A systematic uncertainty related to the knowledge of the trigger efficiencies has been previously studied using large samples of B + → (J/ψ → µ + µ − ) K + and B + → (ψ(2S) → µ + µ − ) K + decays by comparing the ratios of the trigger efficiencies in data and simulation [111]. Based on this comparison, a relative uncertainty of 1.1% is assigned.
Another possible source of uncertainty is the potential disagreement between data and simulation in the estimation of efficiencies due to effects not considered above. This is studied by varying the selection criteria in ranges that lead to changes in the measured signal yields as large as ±20%. For this study, the B + → (ψ(2S) → J/ψη) K + data sample is used. The resulting difference in data-simulation efficiency ratio does not exceed 4.0%, which is conservatively taken as systematic uncertainty.
The systematic uncertainties discussed above affect the ratio of the total efficiencies R ε , and are accounted for in the fits using Gaussian constraints. A different class of systematic uncertainties directly affects the fit itself, namely uncertainties associated with the fit models used to describe the J/ψη and J/ψηK + mass spectra. The systematic uncertainty is accounted for by using fits with alternative models. The alternative resolution models for the C ψ(2S) and C X components include a generalised Student's t-distribution [112] and a sum of two modified Gaussian functions with a power-law tail on each side of the distribution. For wide charmonia and charmonium-like resonances a P-wave relativistic Breit-Wigner function is also tested instead of the S-wave profile. The tail parameters of the F S resolution functions are varied within their uncertainties, as determined from the simulation. For the C NR component, the degree of the polynomial function is varied between zero and two. For the signal component of the fit to the J/ψηK + mass spectrum, the list of alternative models consists of a bifurcated generalised Student's t-distribution [112], an Apollonius function [113], and a sum of two modified Gaussian functions with a power-law tail on each side of the distribution. For the background component, the second-degree positive-definite polynomial function, and the product of an exponential function and a second-degree positive-definite polynomial function, are tested as alternative models. For each alternative model a fit to the J/ψη mass spectrum is performed and the upper limit (UL) on the F X or B X value is determined and conservatively the largest value of the upper limit is taken to account for the systematic uncertainty. For the B + → (ψ 2 (3823) → J/ψη) K + and B + → (ψ(4040) → J/ψη) K + signals, the maximal deviation relative to the baseline fit is taken as uncertainty and added in quadrature to the uncertainty obtained from the fit. For the ψ 2 (3823) state, only 90 (95)% confidence level (CL) upper limits on the natural width of 5.2 (6.6) MeV are known [30]. For this case fits with the natural width varied between 0.2 and 6.6 MeV are performed and the maximal deviation with respect to the default fit, where 1 MeV is assumed, is taken as the corresponding systematic uncertainty.

Results and summary
The upper limits at 90% CL on the ratio of branching fractions F X (m X ) for B + → (X → J/ψη) K + decays via a narrow intermediate X state are set for masses of the hypothetical X particle between 3.7 and 4.7 GeV/c 2 . The upper limits, shown in Fig. 6, are set with the CL s method [114] in which the p-values are calculated based on the asymptotic properties of the profile likelihood ratio [115]. The corresponding upper limits on the product of branching fractions, B X , are calculated in a similar way and shown in Fig. 7. The local statistical significance for the mass values with the weakest upper limits, e.g. m X = 3.952, 4.352 or 4.442 GeV/c 2 , is estimated using Wilks' theorem [116] and is found to be less than three standard deviations.
The asymmetric uncertainty in F ψ 2 (3823) arises from varying the ψ 2 (3823) natural width between 0.2 and 6.6 MeV. The corresponding products of branching fractions are where the second uncertainty is due to the imprecise knowledge of the B + → ψ(2S)K + and ψ(2S) → J/ψη branching fractions [34]. For decays with other intermediate states no signals are seen and the corresponding upper limits are listed in Table 2.
In conclusion, a search for charmonium and charmonium-like exotic states contributing to the J/ψη mass spectrum from B + → J/ψηK + decays is performed, using a data sample corresponding to an integrated luminosity of 9 fb −1 collected with the LHCb detector at 7, 8 and 13 TeV centre-of-mass energies in proton-proton collisions. The B + → (ψ(2S) → J/ψη) K + decay mode is used a normalisation channel. While no narrow resonances are seen, evidence is found for the B + → (ψ 2 (3823) → J/ψη) K + and B + → (ψ(4040) → J/ψη) K + decays, and the corresponding branching fractions and their ratios relative to the normalisation decay mode are measured.