Study of the $\psi_2(3823)$ and $\chi_{c1}(3872)$ states in $B^+ \rightarrow \left( J\psi\pi^+\pi^-\right)K^+$ decays

The decays $B^+\rightarrow J/\psi \pi^+ \pi^- K^+$ are studied using a data set corresponding to an integrated luminosity of 9fb$^{-1}$ collected with the LHCb detector in proton-proton collisions between 2011 and 2018. Precise measurements of the ratios of branching fractions with the intermediate $\psi_2(3823)$, $\chi_{c1}(3872)$ and $\psi(2S)$ states are reported. The decay of $B^+\rightarrow \psi_2(3823)K^+$ with $\psi_2(3823)\rightarrow J\psi\pi^+\pi^-$ is observed for the first time with a significance of 5.1 standard deviations. The mass differences between the $\psi_2(3823)$, $\chi_{c1}(3872)$ and $\psi(2S)$ states are measured to be $$ \begin{array}{rcl} m_{\chi_{c1(3872)}} - m_{\psi_2(3823)}&=&47.50 \pm 0.53 \pm 0.13\,\mathrm{MeV/}c^2\,, \\ m_{\psi_2(3823)} - m_{\psi(2S)}&=&137.98 \pm 0.53 \pm 0.14\,\mathrm{MeV/}c^2\,, \\ m_{\chi_{c1}(3872)} - m_{\psi(2S)}&=&185.49 \pm 0.06 \pm 0.03\,\mathrm{MeV/}c^2\,, \end{array} $$ resulting in the most precise determination of the $\chi_{c1}(3782)$ mass. The width of the $\psi_2(3823)$ state is found to be below 5.2MeV at 90\% confidence level. The Breit-Wigner width of the $\chi_{c1}(3872)$ state is measured to be $$ \Gamma^{\mathrm{BW}}_{\chi_{c1}(3872)} = 0.96^{+0.19}_{-0.18}\pm0.21 \mathrm{MeV},$$ which is inconsistent with zero by 5.5 standard deviations.


Introduction
The observation of a narrow χ c1 (3872) state in the J/ψ π + π − mass spectrum of B + → J/ψ π + π − K + decays by the Belle collaboration in 2003 [1] has led to a renewed interest in the study of hadrons containing heavy quarks. Many new charmonium-like states have since been observed [2]. Some of the new states are unambiguously interpreted as conventional cc states, some are manifestly exotic [3][4][5][6][7][8][9], while for the others a definite interpretation is still missing [10][11][12]. Despite the large amount of experimental data , the nature of the χ c1 (3872) state is still unclear. Several interpretations have been proposed, such as a conventional χ c1 (2P) state [41], a molecular state [42][43][44], a tetraquark [45], a ccg hybrid state [46], a vector glueball [47] or a mixed state [48,49]. Precise measurements of the resonance parameters, namely the mass and the width, are crucial for the correct interpretation of the state. Comparison of the decays of beauty hadrons with final states involving the χ c1 (3872) particle and those involving other charmonium resonances can shed light on the production mechanism, in particular, on the role of D 0 D * 0 rescattering [50].
In this paper, a sample of B + → (X cc → J/ψ π + π − ) K + decays 2 is analysed, where X cc denotes the ψ 2 (3823), χ c1 (3872) or ψ(2S) state and the J/ψ meson is reconstructed in the µ + µ − final state. The study is based on proton-proton (pp) collision data, corresponding to an integrated luminosity of 1, 2, and 6 fb −1 , collected with the LHCb detector at centre-of-mass energies of 7, 8, and 13 TeV, respectively. This data sample allows studies of the properties of the ψ 2 (3823) and χ c1 (3872) states produced in B decay recoiling against a kaon. The presence of the ψ(2S) state in the same sample provides a convenient sample for normalisation and reduction of potential systematic uncertainties. A complementary measurement using inclusive b → (χ c1 (3872) → J/ψ π + π − ) X decays and a data set, corresponding to an integrated luminosity of 1 and 2 fb −1 , collected at the centre-of-mass energies of 7 and 8 TeV, is reported in Ref. [61]. This gives a determination of the resonance parameters for the χ c1 (3872) state with an unprecedented precision, including searches for the poles of the complex Flatté-like amplitude.
classifier based on a decision tree with gradient boosting (BDT) [94].
Muon, pion and kaon candidates are identified by combining information from the RICH, calorimeter and muon detectors [95]. The transverse momentum of muon (hadron) candiates is required to be larger than 550 (220) MeV/c. To allow for efficient particle identification, kaons and pions are required to have a momentum between 3.2 and 150 GeV/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.0 and 3.2 GeV/c 2 .
To form the B + candidates, the selected J/ψ candidates are combined with a pair of oppositely charged pions and a positively charged kaon. 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 particle under consideration. To improve the mass resolution for the B + candidates, a kinematic fit [96] is performed. This fit constrains the mass of the µ + µ − pair to the known mass of the J/ψ meson [2] and constraints the B + candidate to originate from its associated PV. In addition, the measured decay time of the B + candidate, calculated with respect to the associated PV, is required to be greater than 75 µm/c. This requirement suppresses background from particles originating from the PV.
A BDT is used to further suppress the combinatorial background. It is trained using a simulated sample of B + → (ψ 2 (3823) → J/ψ π + π − ) K + decays as the signal. For the background, a sample of J/ψ π + π + K − combinations with same-sign pions in data, passing the preselection criteria and having the mass in the range between 5.20 and 5.35 GeV/c 2 , is used. The k-fold cross-validation technique [97] with k = 13 is used to avoid introducing a bias in the BDT evaluation. The BDT is trained on variables related to the reconstruction quality, decay kinematics, decay time of B + candidate and the quality of the kinematic fit. The requirement on the BDT output is chosen to maximize ε/(α/2 + √ B) [98], where ε is the signal efficiency for the B + → ψ 2 (3823)K + decays obtained from simulation; α = 5 is the target signal significance in units of standard deviations; B is the expected background yield within narrow mass windows centred at the known B + and ψ 2 (3823) masses [2]. The mass distribution of selected B + → J/ψ π + π − K + candidates is shown in Fig. 1. The data are fit with a sum of a modified Gaussian function with power-law tails on both sides [99,100] and a linear polynomial combinatorial background component. The B + signal yield is (547.8 ± 0.8) × 10 3 candidates.
For each of the three regions the J/ψ π + π − K + mass is restricted to 5.20 ≤ m J/ψ π + π − K + < 5.35 GeV/c 2 . To improve the resolution on the J/ψ π + π − mass and to eliminate a small correlation between m J/ψ π + π − K + and m J/ψ π + π − variables, the m J/ψ π + π − variable is computed using a kinematic fit [96] that constrains the mass of the B + candidate to its known value [2]. In each region, the fit function is defined as a sum of four components: 1. signal B + → X cc K + decays parameterised as a product of the B + and X cc signal templates described in detail in the next paragraph; 2. contribution from the decays B + → (J/ψ π + π − ) NR K + with no narrow intermediate X cc state, parameterised as a product of the B + signal template and a linear function of m J/ψ π + π − ; 3. random combinations of X cc and K + candidates, parameterised as a product of the X cc signal template and a linear function of m J/ψ π + π − K + ; 4. random J/ψ π + π − K + combinations, described below.
The templates for the B + signals are described by a modified Gaussian function with power-law tails on both sides of the distribution [99,100]. The tail parameters are fixed to the values obtained from simulation. The narrow X cc signal templates are parameterised with S-wave relativistic Breit-Wigner functions convolved with the mass resolution. Due to the proximity of the χ c1 (3872) state to the D 0 D * 0 threshold, modelling this component as a Breit-Wigner function may not be adequate [101][102][103][104][105]. However, the analysis from Ref. [61] demonstrates that a good description of data is obtained with a Breit-Wigner lineshape when the mass resolution is included. The mass resolution is described by a symmetric modified Gaussian function with power-law tails on both sides of the distribution, with the parameters fixed to the values from simulation. In the template for the B + signal, the peak-position parameter is shared between all three decays and allowed to vary in the fit. The mass resolutions used in the B + and X cc signal templates are fixed to the values determined from simulation, but are corrected by common scale factors, f B + and f X cc , to account for a small discrepancy in the mass resolution between data and simulation. The masses of the X cc signal templates, as well as the Breit-Wigner widths for the ψ 2 (3823) and χ c1 (3872) states, are free fit parameters, while the width in the template for the ψ(2S) signal is fixed to its known value [2]. The combinatorial-background component is modelled with a smooth two-dimensional function where E(m J/ψ π + π − K + ) is an exponential function, P 3,4 (m J/ψ π + π − ) is a three-body phase-space function [106], and P 2D is a two-dimensional positive bilinear function, which accounts for small non-factorizable effects. For the considered fit ranges P 3,4 (m J/ψ π + π − ) is close to a constant. The J/ψ π + π − K + and J/ψ π + π − mass distributions together with projections of the simultaneous unbinned maximum-likelihood fit are shown in Fig. 2. Signal yields N B + →X cc K + , calculated mass differences δm X cc ≡ m X cc − m ψ(2S) , Breit-Wigner widths Γ X cc and resolution scale factors are listed in Table 1. The fit model is tested using pseudoexperiments and no bias is found in the results and their associated uncertainties. The masses of B + and ψ(2S) mesons are found to be compatible with their known values [2]. The fit component corresponding to the B + → (J/ψ π + π − ) NR K + is found to be negligible for the ψ(2S) region, dominant for the ψ 2 (3823) region and small for the χ c1 (3872) region. The fit component corresponding to the random X cc K + combinations is negligible for all fit regions. The statistical significance of the observed B + → (ψ 2 (3823) → J/ψ π + π − ) K + signal over the background-only hypothesis is estimated to be 5.1 standard deviations using Wilks' theorem [107]. The significance is confirmed by simulating a large number of pseudoexperiments according to the background distribution observed in data.
The likelihood profiles for the Breit-Wigner widths of ψ 2 (3823) and χ c1 (3872) states are presented in Fig. 3. From these profiles the Breit-Wigner width of the χ c1 (3872) state is found to be inconsistent with zero by 5.5 standard deviations, while for the ψ 2 (3823) state the width is consistent with zero.
shown as points with error bars. A fit, described in the text, is overlaid.

Ratios of branching fractions
Ratios of branching fractions, R X Y , are defined as  where X, Y stand for either the ψ 2 (3823), χ c1 (3872) or ψ(2S) states. They are estimated as where N is the signal yield reported in Table 1 and ε denotes the efficiency of the corresponding decay. The efficiency is defined as the product of geometric acceptance, reconstruction, selection, hadron identification and trigger efficiencies, where each subsequent efficiency is defined with respect to the previous one. All of the contributions, except that of the hadron-identification efficiency, are determined using simulated samples. The hadron-identification efficiency is determined using large calibration samples of for kaons and pions [69,108]. The ratios of the efficiencies are determined to be where the uncertainty reflects the limited size of the simulated samples. Other sources of systematic uncertainty are discussed in the following section. The ratios of the efficiencies differ from unity mostly due to the different pion momentum spectra in the different X cc → J/ψ π + π − decays.

Systematic uncertainty
Due to the similar decay topologies, systematic uncertainties largely cancel in the ratios R X Y . The remaining contributions are listed in Table 2 and are discussed below. The systematic uncertainty related to the signal and background shapes is investigated using alternative parameterisations. A generalized Student's t-distribution [109], an Apollonios function [110] and a modified Novosibirsk function [111] are used as alternative models for the B + signal template. For the X cc signal template, alternative parameterisations of the mass resolution, namely a symmetric variant of an Apollonios function [110], a Student's t-distribution and a sum of two Gaussian functions sharing the same mean are considered. In addition, P-wave and D-wave relativistic Breit-Wigner functions are used as alternative ψ 2 (3823) signal templates, and the Blatt-Weisskopf barrier factors [112] are varied between 1.5 and 5 GeV −1 . The width of the ψ(2S) state, fixed in the fit, is varied between 270 and 302 keV [2]. The maximal deviations in the ratios R X Y with respect to the baseline fit model are taken as systematic uncertainties for each of the systematic signal model sources. For the systematic uncertainty related to the modelling of the smooth polynomial functions, pseudoexperiments with about 10 7 simulated events (approximately 100 times large than data sample are generated with the baseline fit model and fitted with alternative background models. In this study the degree of the polynomial functions is varied from the first to the second order, separately for each fit component and each channel. In each case the ratio R X Y is computed and the maximal difference with respect to the baseline fit model is taken as a corresponding systematic uncertainty. Since the decay model for ψ 2 (3823) → J/ψ π + π − is unknown, a phase-space model is used in simulation. To probe the associated systematic uncertainty the model discussed in Ref. [59] is used. This model accounts for the quantum-chromodynamics multipole expansion [113], as well as the effective description of the coupled-channel effects via hadronic-loop mechanism [114] with the interference phase Φ as a free parameter. The π + π − mass spectrum and the angular distributions in the decay strongly depend on the phase Φ, however, the efficiency for the B + → (ψ 2 (3823) → J/ψ π + π − ) K + decays is found to be stable. It varies within 0.2% with respect to the efficiency computed for the phase-space model when the unknown phase Φ varies in the range −π ≤ Φ < π.
An additional uncertainty arises from differences between the data and simulation, in particular differences in the reconstruction efficiency of charged-particle tracks.
The track-finding efficiencies obtained from the simulation samples are corrected using data-driven techniques [82]. The uncertainties related to the correction factors, together with the uncertainty in the hadron-identification efficiency due to the finite size of the calibration samples [69,108], are propagated to the ratio of total efficiencies using pseudoexperiments.
The systematic uncertainty related to the trigger efficiency is estimated using large samples of the B + → J/ψ K + and B + → ψ(2S)K + decays by comparing the ratios of trigger efficiencies in data and simulation [83]. The imperfect data description by the simulation due to remaining effects is studied by varying the BDT selection criteria in ranges that lead to ±20% change in the measured efficiency. The resulting variations in the efficiency ratios do not exceed 1%, which is taken as a corresponding systematic uncertainty. The last systematic uncertainty considered for the ratio R X Y is due to the finite size of the simulated samples.
For each choice of the fit model, the statistical significance of the observed B + → (ψ 2 (3823) → J/ψ π + π − ) K + signal is calculated from fit to data using Wilks' theorem. The smallest significance found is 5.1 standard deviations, numerically close to the value obtained from the baseline fit model.
The systematic uncertainties on the mass differences between the ψ 2 (3823), χ c1 (3872) and ψ(2S) states are summarized in Table 3. An important source of systematic uncertainty is due to the signal and background shapes. Different parameterisations of the signal templates and non-signal components, described above, are used as the alternative fit models. The maximal deviation in the mass differences with respect to the baseline results is assigned as the corresponding systematic uncertainty. The uncertainty in the momentum-scale calibration, important for mass measurements, e.g. Refs. [51,67,68,84, 89, 92, 93, 115-124], largely cancels for the mass differences. The associated systematic uncertainty is evaluated by varying the momentum scale within its known uncertainty [68] and repeating the fit. The J/ψ π + π − mass is computed constraining the mass of the B + candidate to the known value, m B + = 5279.25 ± 0.26 MeV/c 2 [2]. The uncertainty on the B + meson mass is propagated to the measured mass differences.
The main source of systematic uncertainty for the Breit-Wigner widths Γ ψ 2 (3823) and The systematic uncertainty due to the mismodelling of the experimental resolution in simulation is accounted for with the resolution scale factors f B + and f X cc and therefore is included as a part of the statistical uncertainty. A small dependency of the scale factor f X cc on the dipion momentum for the ψ(2S) → J/ψ π + π − decay is reported in Ref. [61]. Such effect causes a bias in the effective scale factor for different decays due to slightly different dipion spectra. Such bias is found to be negligible with respect to the statistical uncertainty for the factor f X cc . The analysis is carried out by neglecting any interference effects between the X cc resonances and other components. Such an assumption can bias the measurement of the mass and width-parameters associated to the X cc states. To account for such interference effects a full amplitude analysis is required, which is beyond the scope of this study. However, to estimate the possible effect of this assumption on the χ c1 (3872) mass and width-parameters, the background-subtracted J/ψ π + π − mass distribution in the χ c1 (3872) region is studied with the sPlot technique used for background subtraction [125] using the J/ψ π + π − K + mass as the discriminative variable. The distribution is fit with a model that accounts for the signal, coherent and incoherent backgrounds where A BW (m) is a Breit-Wigner amplitude, convolved with the mass resolution function R, and N stands for a normalisation constant. The coherent and incoherent background components b c (m) and b 2 i (m) are parameterised with polynomial functions. The relative interference phase δ(m) is taken to be constant for the narrow 3.85 ≤ m J/ψ π + π − < 3.90 GeV/c 2 region, δ(m) ≡ δ 0 . An equally good description of data is achieved for totally incoherent (b c (m) ≡ 0) and coherent (b 2 i (m) ≡ 0) background hypotheses, as well as for any intermediate scenarios with the phase δ 0 close to π 2 . The latter reflects a high symmetry of the observed χ c1 (3872) lineshape. For all scenarios, variations of the mass and width parameters are limited to 50 keV/c 2 and 150 keV, respectively.

Results and summary
The decay of B + → (ψ 2 (3823) → J/ψ π + π − ) K + is observed for the first time with a significance of 5.1 standard deviations.
The value of the Breit-Wigner width Γ χ c1 (3872) agrees well with the value from the analysis of a large sample of χ c1 (3872) → J/ψ π + π − decays from the inclusive decays of beauty hadrons [61]. Using the known value of the ψ(2S) mass [2], the Breit−Wigner masses for the ψ 2 (3823) and χ c1 (3872) states are computed to be m ψ 2 (3823) = 3824.08 ± 0.53 ± 0.14 ± 0.01 MeV/c 2 , m χ c1 (3872) = 3871.59 ± 0.06 ± 0.03 ± 0.01 MeV/c 2 , where the last uncertainty is due to the knowledge of the ψ(2S) mass. These are the most precise measurements of these masses. The mass difference between χ c1 (3872) and ψ(2S) states is more precise than the average reported in Ref. [2]. It also agrees well with the measurement from Ref. [61]. Taking into account a partial overlap of the data sets and correlated part of systematic uncertainty, the LHCb average mass difference and the mass of the χ c1 (3872) state are where the second uncertainty is due to the knowledge of the ψ(2S) mass. The difference between the m χ c1 (3872) mass, determined from the Breit-Wigner fit, and the D 0 D * 0 threshold δE ≡ ( m D 0 + m D * 0 ) c 2 − m χ c1 (3872) c 2 is computed to be δE = 0.12 ± 0.13 MeV , δE| LHCb = 0.07 ± 0.12 MeV , where the first value corresponds to the measurement performed in this analysis, while the second one is an average with results from Ref. [61]. A value of 3871.70 ± 0.11 MeV/c 2 is taken for the threshold m D 0 + m D * 0 , calculated from Ref. [2,61], accounting for the correlation due to the knowledge of the charged and neutral kaon masses between the measurements. The uncertainty on δE is now dominated by the knowledge of kaon masses. These are the most precise measurements of the χ c1 (3872) mass and δE parameter. [25] BaBar collaboration, P.  [89] LHCb collaboration, R. Aaij et al., Observation of Λ 0 b → ψ(2S)pK − and Λ 0 b → J/ψ π + π − pK − decays and a measurement of the Λ 0 b baryon mass, JHEP 05 (2016) 132, arXiv:1603.06961.