Observation of sizeable $\omega$ contribution to $\chi_{c1}(3872)\to\pi^+\pi^-J/\psi$ decays

Resonant structures in the dipion mass spectrum from $\chi_{c1}(3872)\to\pi^+\pi^- J/\psi$ decays, produced via $B^+\to K^+\chi_{c1}(3872)$ decays, are analyzed using proton-proton collision data collected by the LHCb experiment, corresponding to an integrated luminosity of 9 $fb^{-1}$. A sizeable contribution from the isospin conserving $\chi_{c1}(3872)\to\omega J/\psi$ decay is established for the first time, $(21.4\pm2.3\pm2.0)\%$, with a significance of more than $7.1\sigma$. The amplitude of isospin violating decay, $\chi_{c1}(3872)\to\rho^0 J/\psi$, relative to isospin conserving decay, $\chi_{c1}(3872)\to\omega J/\psi$, is properly determined, and it is a factor of six larger than expected for a pure charmonium state.

(background) shape modelled by a double-sided Crystal Ball [17] (quadratic) function, yielding 6788 ± 117 χ c1 (3872) → π + π − J/ψ decays with a χ c1 (3872) mass resolution of σ m = 2.66 ± 0.09 MeV (natural units are used throughout). The signal purity is 77% within ±2 σ m around the χ c1 (3872) mass. The dominant source of background is from B + decays to J/ψ meson and excited kaons (K * + ), which decay to K + π + π − . The dipion mass distribution is obtained by two-dimensional unbinned fits of the χ c1 (3872) signal yields to the (m π + π − J/ψ , m π + π − ) data in m π + π − intervals. The signal shape in m π + π − J/ψ is fixed to the global fit result of Fig. 1, while the background shape may vary in each m π + π − interval. The signal and the background m π + π − J/ψ shapes are then multiplied by a two-body phase-space factor, i.e. the J/ψ momentum in the χ c1 (3872) rest frame (p J/ψ ), which depends on both m π + π − J/ψ and m π + π − . The resultant shapes are normalized to unity using their integral over the fitted phase-space range. The obtained data points are shown in Fig. 2. Signal simulation is used to obtain the m π + π − dependence of the dipion mass resolution and of the reconstruction efficiency. In the simulation, pp collisions are generated using Pythia [18] with a specific LHCb configuration [19]. Decays of unstable particles are described by EvtGen [20], in which final-state radiation is generated using Photos [21]. In particular, B + → K + χ c1 (3872), χ c1 (3872) → J/ψρ 0 , ρ 0 → π + π − , J/ψ → µ + µ − signal decays are simulated using the helicity model where the χ c1 (3872) decays via an S-wave, which describes the angular distributions in the data well [5]. When integrating over m π + π − , the generated m π + π − distribution is weighted to reproduce the observed distribution. The interaction of the generated particles with the detector, and its response, are implemented using the Geant4 toolkit [22] as described in Ref. [23]. The underlying pp interaction is reused multiple times, with an independently generated signal decay for each event [24]. The transverse-momentum distribution of the signal B + is weighted to match the data. The relative signal reconstruction efficiency obtained from the simulation is well described by a quadratic function, where m π + π − and σ(m π + π − ) are in MeV. The simulation underestimates the χ c1 (3872) and B + mass resolution by 6% and 14%, respectively. Therefore, σ(m π + π − ) is scaled up by f m = 1.06, which is varied from 1.00 to 1.14 to assess a systematic uncertainty. All theoretical probability density functions fit to the data, P(m π + π − ), are multiplied by the relative efficiency function, ϵ(m π + π − ) and convolved with a Gaussian distribution characterized by a root-mean-square value of f m σ(m π + π − ). The matrix element, M, describing the three-body decay, χ c1 (3872) → π + π − J/ψ, is related to P(m π + π − ) via a scaling (S) and the phase-space factors, P(m π + π − ) = S p J/ψ p π |M| 2 , where p π is the pion momentum in the rest frame of the π + π − system. The scaling factor S is a nuisance parameter.
Following the previous analyses [3,13], the ρ 0 component is first parameterized as a Breit-Wigner (BW) amplitude, ] is the Blatt-Weisskopf barrier factor for the P -wave decay of a vector particle to π + π − , and contains an effective hadron-size parameter R. With R adjusted to 1.45 GeV −1 , a complex phase of the BW amplitude varies in the fitted m π + π − range within one degree of the isovector π + π − P -wave parametrization extracted from the scattering data by the phenomenological analysis of Ref. [25]. A similar value of 1.5 GeV −1 was used in the previous analyses [3,13]. The χ 2 fit, with S as the only fit parameter, fails to describe the data as shown in Fig. 2, with a χ 2 value per number of degrees of freedom (χ 2 /NDoF) equal to 366.6/34. A large disagreement between the data and the χ c1 (3872) → ρ 0 J/ψ amplitude at high dipion masses was missed in the previous comparisons with the χ c1 (3872) → ρ 0 J/ψ simulation (see Fig. S2 in the supplementary material of Ref. [5]; see also Refs. [26,27]), because the EvtGen [20] model does not correctly simulate the impact of the phase-space factors (here p J/ψ p π ) on resonant shapes. As a consequence, the large ρ 0 − ω interference in the data (see below) was mistakenly interpreted as a part of the ρ 0 resonance itself.
Adding an ω contribution via the BW sum model with M = BW(s|m ρ , Γ ρ ) + a exp(i ϕ) BW(s|m ω , Γ ω ), where a and ϕ are the amplitude and phase of the ω component relative to the ρ 0 term [3,13], improves the fit quality substantially, χ 2 /NDoF = 102.9/33, yet not enough to be acceptable. Summing BW amplitudes results in matrix elements that are not unitary, violating first principles of scattering theory. A two-channel K-matrix (K) model [28], coupling the π + π − and π + π − π 0 channels via an ω contribution, resolves this issue, where g are the coupling constants described later. The g 2 ρ→3π coupling is 4-5 orders of magnitude smaller than g 2 ρ→2π [10,29,30] and has been neglected here. The T -matrix is obtained fromT = [1 − i K ϱ] −1 K, where the phase-space matrix ϱ is diagonal, ϱ = diag(ϱ 2π (s), ϱ 3π (s)), and is detailed in the supplemental material. The decay amplitude is given by whereÂ 2π ≡ α 2πT2π,2π + α 3πT2π,3π , and the elements of the production Q-vector (α 2π , α 3π ) are real [31]. The coupling constants are fully determined from other experiments [10], Numerically, g 2 ω→2π /g 2 ρ→2π ≈ 0.0009, while g ω→2π g ω→3π /g 2 ρ→2π ≈ 0.01. Thus, the diagonal ω coupling to two pions can be neglected in comparison with the off-diagonal coupling. In this approximation, equivalent to the full formalism given the precision of this analysis, which is the ρ 0 Breit-Wigner amplitude (Eq. 1) and α 2π is the ρ 0 production factor in the χ c1 (3872) decay. A mild dependence of α 2π on s is possible, since α 2π (s) must be analytic and cannot change much within the resonance widths (see e.g. Ref. [32]). Using Chebyshev polynomials (C n ) to express the dependence: α 2π (s) = n=N n=0 P n C n (ŝ), whereŝ ≡ 2 (s − s min )/(s max − s min ) − 1, s min = (380 MeV) 2 and s max = (775 MeV) 2 . The polynomial coefficients P n are determined from a fit to the data, except for P 0 , which is set to unity as the normalization choice. For α 2π (s) to have the expected theoretical behavior, the series must converge, |P n+1 | < |P n | and N should be kept small.
The ω contribution enters via the element, This term approaches zero at the bare ρ 0 pole mass, which complicates probing the ω contribution, m ω = 782.66 ± 0.13 MeV and Γ ω = 8.68 ± 0.13 MeV [10]. This zero is an artifact of the K-matrix approach, rather than an expectation from scattering theory. To remove it and restore a more physical behavior of the ω term, is set. The constant term m ω 2 − m ρ 2 is introduced above to make a ω express the ratio of ω/ρ amplitudes at the ω pole. The value of a ω is determined by the fit to the data. Using Eqs. 7, 8 and 9, with α 2π constant (N = 0), is the common method to describe ρ 0 − ω interference in analyses of the π + π − system [33]. The exact K-matrix formulae are used here. In view of Eq. 7, models with α 3π = 0 are interpreted as containing a ρ 0 component only. The following integrals are calculated to quantify a relative rate of the ω contribution: I tot = P(m π + π − )dm π + π − , I ρ = P(m π + π − | α 3π = 0)dm π + π − and I ω = P(m π + π − | α 2π = 0)dm π + π − , where P(m π + π − ) is neither convolved with the mass resolution, nor multiplied by the efficiency function. The integration is performed over the full phase space. To quantify the overall impact of ω on the total rate, including ρ 0 − ω interference effects, the ratio R all ω ≡ 1 − I ρ /I tot is defined. The traditional ω fit fraction is given by R 0 ω ≡ I ω /I tot . Finally, to probe the ratio of the χ c1 (3872) isospin conserving to isospin violating couplings, the ratio R 0 ω/ρ ≡ I ω /I ρ is introduced. Using the B + → K + χ c1 (3872) simulation, the extraction of the dipion mass distribution and subsequent fit with the K-matrix model, with the efficiency and the mass resolution corrections included, are verified to a precision an order of magnitude better than the total systematic uncertainties reported below.
A number of analysis variations have been performed to evaluate systematic uncertainties on the ω fraction, as summarized in Table 1. In the dipion mass extraction from the data, Gaussian functions are used for the χ c1 (3872) shape in m π + π − J/ψ data, rather than Crystal Ball functions. A cubic polynomial replaces the quadratic function in the relative efficiency parametrization.
Data-driven corrections to the simulation of the hadron identification are implemented, as a cross-check. To check further the relative efficiency simulation and the background subtraction, a tighter data selection is implemented with a Boosted Decision Tree (BDT) classifier [35], using inputs from hadron identification variables, the B + vertex quality, the PV impact parameters of the final-state particles and the B + candidate, the hadron transverse-momenta and the B + flight distance. The tighter selection reduces the combinatoric background under the B + peak in the m K + π + π − J/ψ distribution from 9.4% to 3.1%, resulting in an overall background reduction under the χ c1 (3872) peak in the m π + π − J/ψ distribution from 23% to 17%. At the same time, the χ c1 (3872) signal yield is reduced by only 0.4%.
To evaluate uncertainties due to the mass resolution used in the fit, the scaling factor f m is varied and the fit range is reduced to exclude the last interval of the m π + π − distribution (775-780 MeV), since this falls beyond the phase-space limit, and therefore its content is very sensitive to σ(m π + π − ).
To check for interference effects between B + → K + χ c1 (3872) and B + → K * + J/ψ decays, the data are split into two subsamples depending on the sign of cos θ X , where θ X is the χ c1 (3872) helicity angle, the angle between the K + and J/ψ momenta in the χ c1 (3872) rest frame. The composition of K * + resonances is different in each subsample.
As a variation of the production model, non-resonant (NR) terms are added to the production vector, Without χ c1 (3872) → 3πJ/ψ data in the fit, the NR production parameter f 3π cannot be probed, thus it is set to zero. A good-quality fit is obtained with constant α 2π , which gives f 2π = (−9.7 ± 1.6) × 10 −7 .
The default fit assumes an S-wave χ c1 (3872) decay. When a D-wave component is allowed in the fit, multiplying |M| 2 by 1+A D 2 B 2 (p J/ψ )/B 2 (p ρ J/ψ ), where p ρ J/ψ ≡ p J/ψ (m ρ ) and B 2 (p) ≡ p 4 /(9 + 3 (R p) 2 + (R p) 4 ), the A D parameter is consistent with zero, 0.13 ± 0.41. Tuning A D to 0.176 produces a 4% D-wave fraction, equal to the upper limit from studies of the angular correlations [5]. Since the hadron size parameter R is tuned to the scattering data, which automatically includes any ρ 0 excitations, there is no strong motivation to include the ρ(1450) pole [10] (ρ ′ ) in the K-matrix model. When included, the fit quality is slightly reduced.
The size parameter R used in B 1 (p π ) of Eq. 3 could be related to the χ c1 (3872) size, rather than the ρ 0 or ω sizes. Thus, fit variations are tried in which it has an independent value, R prod , varied from zero to 30 GeV −1 . The R value itself is varied within the range 1.3-1.6 GeV −1 . An alternative model of the ρ 0 shape, that does not include an R-dependent Blatt-Weisskopf form factor, is provided by the Gounaris-Sakurai (GS) formula [36], and was utilized by the BaBar collaboration to describe high statistics e + e − → π + π − (γ) data [37]. As discussed in the supplemental material, when applying this prescription with ρ 0 , ρ ′ and ω contributions, an excellent fit to the data is obtained, χ 2 /NDoF = 24.8/32, p-value = 81%, matching the fit quality of the default fit.
Total systematic uncertainties, given at the bottom of Table 1, are set to cover the maximal deviation from the default fit results, and are comparable to the statistical uncertainties. The lowest ω significance is 7.1σ, excluding the subsamples and P 2 ̸ = 0 fit as discussed above. This is a more significant observation of χ c1 (3872) → ωJ/ψ decays than achieved using the dominant ω → π + π − π 0 decay channel.
In addition to the resonant coupling constants, the limited phase space in χ c1 (3872) → π + π − J/ψ decays has a large impact on the R 0 ω/ρ value capturing a much smaller fraction of the ω resonance than for the ρ 0 resonance. To probe the resonant coupling constants, the phase space can be artificially extended to contain both resonance peaks by setting m χ c1 (3872) = 4000 MeV, as illustrated in the supplemental material. Integrating the default P(m π + π − ) model in the extended phase space, R 0 ω/ρ ′ = 0.18 ± 0.05 is obtained and then used to deduce the ratio of the isospin violating to isospin conserving χ c1 (3872) couplings, This value is an order of magnitude larger than expected for pure cc states, as exemplified by the corresponding value for the ψ(2S) state [10], where p η and p π 0 are the breakup momenta [38]. Therefore, the χ c1 (3872) state cannot be a pure charmonium state. Large isospin violation is naturally expected in models in which the χ c1 (3872) state has a significant DD * component, either as constituents (i.e. in the "molecular model") or generated dynamically in the decay [39][40][41][42][43][44][45][46][47][48]. The proximity of the χ c1 (3872) mass to the D 0 D * 0 threshold, enhances such contributions over D + D * − combinations, which are 8 MeV heavier. It has also been suggested in compact tetraquark models that two neutral states could be degenerate and mix, giving rise to large isospin violation in χ c1 (3872) decays [49][50][51].
In summary, the ρ 0 and ω contributions to χ c1 (3872) → π + π − J/ψ decays are resolved for the first time using a much larger data sample than previously available. Through ρ 0 − ω interference, the ω contribution accounts for (21.4 ± 2.3 ± 2.0)% of the total rate, or equivalently, B(χ c1 (3872) → ρ 0 J/ψ)/B(χ c1 (3872) → π + π − J/ψ) = (78.6 ± 2.3 ± 2.0)%. Excluding interference effects, the ω contribution, (1.9 ± 0.4 ± 0.3)%, is found to be consistent with, but more precise than the previous χ c1 (3872) → ωJ/ψ measurements utilizing ω → π + π − π 0 decays [9,11,12]. The isospin violating ρ 0 contribution, quantified for the first time with proper subtraction of the ω contribution, is an order of magnitude too large for χ c1 (3872) to be a pure charmonium state. This analysis also serves as an excellent illustration for why a simple sum of ρ 0 and ω Breit-Wigner amplitudes should not be used to describe interfering resonances. Such a model fails to describe the data, while the K-matrix approach, which respects unitarity, provides an excellent description of the data and gives numerical results that are consistent with the previous χ c1 (3872) → ωJ/ψ measurements utilizing ω → π + π − π 0 decays. Observation of sizeable ω contribution to χ c1 (3872) → π + π − J/ψ decays Supplemental material 1 K-matrix phase-space matrix elements A notation in which the Blatt-Weisskopf barrier factors (B l ) are integrated with the phase-space matrix is used in this Letter, (S1) Since the ω meson decays to three pions via a P -wave ρπ decay, 2 Gounaris-Sakurai resonant m π + π − shape An alternative model to the R-dependent Blatt-Weisskopf form factor in the Breit-Wigner amplitude is provided by the Gounaris-Sakurai formula [36], where, , and h ′ (s) is the derivative of h(s), with respect to s, calculated numerically. To include the P -wave momentum barrier in the ρ 0 decay, the matrix element is set to M = p π (s)/p π (m ρ ) BW GS ρ (s, m ρ , Γ ρ ). A fit of this ρ 0 shape to the data is only slightly Gounaris-Sakurai model Figure S1: Distribution of m π + π − in χ c1 (3872) → π + π − J/ψ decays, fit with the Gounaris-Sakurai model. The dashed line represents the ρ ′ contribution multiplied by a factor of 10.
better than the fit of the Breit-Wigner shape given by Eq. 1 (χ 2 /NDoF = 290.0/34, p-value = 2 × 10 −42 ). A good-quality fit to the data is achieved following the prescription used by the BaBar collaboration to describe a large e + e − → π + π − (γ) data sample [37], Following this work, a simple Breit-Wigner amplitude for the ω meson is used, BW ω (s, m ω , Γ ω ) = m 2 ω /(m 2 ω − s − im ω Γ ω ), and the ρ 0 , ω and ρ ′ masses and widths are taken from Table VI of Ref. [37]. The term including the ρ 0 and ω resonances (Eq. S9) is equivalent to Eqs. 8-9, and thus originates from the coupled-channel approach. No complex phase is expected in this approach, thus ϕ ω is set to zero. Using the small phase, consistent with zero, obtained by the BaBar collaboration, gives almost identical results. Since, the ρ ′ contribution is not added via the K-matrix, its phase can be different from zero and is fixed to the central value of the BaBar result, ϕ ρ ′ = 3.76 ± 0.10 rad. The m π + π − distribution, fit with this model, is shown in Fig. S1. The fit quality is excellent (χ 2 /NDoF = 24.8/32, p-value = 0.81), matching the fit quality of the default fit. The ρ ′ significance is 3.1σ using Wilks' theorem [34]. Its production parameter relative to the ρ 0 meson, A GS ρ ′ = 0.302 ± 0.099, is consistent within the large uncertainty with the value obtained by the fit to the BaBar data, 0.158 ± 0.018 [37]. The ω significance is 7.8σ, and its production parameter, A GS ω = 0.0171 ± 0.0024, is an order of magnitude larger than that obtained by the BaBar collaboration, 0.001644 ± 0.00061 [37]. This is not surprising, since in e + e − collisions, the ρ 0 meson is produced via electromagnetic interactions, which do not follow isospin symmetry, while in χ c1 (3872) decays, the ρ 0 state is produced via strong interactions, which suppress isospin violating decays. The ω fractional contributions are consistent with the default model, R all ω = 0.221 ± 0.024, R 0 ω = 0.021 ± 0.005 and R 0 ω/ρ = 0.028 ± 0.007. Here, R all ω = 1 − R ρ+ρ ′ , where R ρ+ρ ′ = 0.780 ± 0.023 is the  Figure S2: Amplitude model for χ c1 (3872) → π + π − J/ψ decays, obtained by the default fit to the LHCb data (Fig. 3), with the phase-space integration limit extended by setting the χ c1 (3872) mass to 4000 MeV. The actual phase-space limit, imposed by the measured χ c1 (3872) mass, is indicated by the vertical dashed line. Neither mass resolution nor detector efficiency are applied here.