Evidence for an η c ( 1 S ) π - resonance in B 0 → η c ( 1 S ) K + π - decays.

A Dalitz plot analysis of B 0 → η c ( 1 S ) K + π - decays is performed using data samples of pp collisions collected with the LHCb detector at centre-of-mass energies of s = 7 , 8 and 13 Te V , corresponding to a total integrated luminosity of 4.7 fb - 1 . A satisfactory description of the data is obtained when including a contribution representing an exotic η c ( 1 S ) π - resonant state. The significance of this exotic resonance is more than three standard deviations, while its mass and width are 4096 ± 20 - 22 + 18 Me V and 152 ± 58 - 35 + 60 Me V , respectively. The spin-parity assignments J P = 0 + and J P = 1 - are both consistent with the data. In addition, the first measurement of the B 0 → η c ( 1 S ) K + π - branching fraction is performed and gives B ( B 0 → η c ( 1 S ) K + π - ) = ( 5.73 ± 0.24 ± 0.13 ± 0.66 ) × 10 - 4 , where the first uncertainty is statistical, the second systematic, and the third is due to limited knowledge of external branching fractions.

Introduction

Since the discovery of the X(3872) state in 2003 [1], several exotic hadron candidates have been observed, as reported in recent reviews [2-7].1 The decay modes of these states indicate that they must contain a heavy quark–antiquark pair in their internal structure; however, they cannot easily be accommodated as an unassigned charmonium or bottomonium state due to either their mass, decay properties or electric charge, which are inconsistent with those of pure charmonium or bottomonium states. Different interpretations have been proposed about their nature [2-4], including their quark composition and binding mechanisms. In order to improve the understanding of these hadrons, it is important to search for new exotic candidates, along with new production mechanisms and decay modes of already observed unconventional states.

The state, discovered by the BESIII collaboration in the final state [9], and confirmed by the Belle  [10] and CLEO  [11] collaborations, can be interpreted as a hadrocharmonium state, where the compact heavy quark–antiquark pair interacts with the surrounding light quark mesonic excitation by a QCD analogue of the van der Waals force [12]. This interpretation of the state predicts an as-yet-unobserved charged charmonium-like state with a mass of approximately whose dominant decay mode is to the system.2 Alternatively, states like the meson could be interpreted as analogues of quarkonium hybrids, where the excitation of the gluon field (the valence gluon) is replaced by an isospin-1 excitation of the gluon and light-quark fields [13]. This interpretation, which is based on lattice QCD, predicts different multiplets of charmonium tetraquarks, comprising states with quantum numbers allowing the decay into the system. The system carries isospin , G-parity , spin and parity , where L is the orbital angular momentum between the and the mesons. Lattice QCD calculations [14, 15] predict the mass and quantum numbers of these states, comprising a state of mass , a state of mass , and a state of mass . The resonance, discovered by the Belle collaboration [16] and confirmed by LHCb  [17, 18], could also fit into this scenario. Another prediction of a possible exotic candidate decaying to the system is provided by the diquark model [19], where quarks and diquarks are the fundamental units to build a rich spectrum of hadrons, including the exotic states observed thus far. The diquark model predicts a candidate below the open-charm threshold that could decay into the final state. Therefore, the discovery of a charged charmonium-like meson in the system would provide important input towards understanding the nature of exotic hadrons.

In this article, the decay is studied for the first time, with the meson reconstructed using the decay mode. The decay is expected to proceed through intermediate states, where refers to any neutral kaon resonance, following the diagram shown in Fig. 1a. If the decay also proceeds through exotic resonances in the system, denoted by states in the following, a diagram like that shown in Fig. 1b would contribute. The decay involves only pseudoscalar mesons, hence it is fully described by two independent kinematic quantities. Therefore, the Dalitz plot (DP) analysis technique [20] can be used to completely characterise the decay.

Fig. 1

Feynman diagrams for a and b decay sequences

The data sample used corresponds to an integrated luminosity of of collision data collected with the LHCb detector at centre-of-mass energies of and in 2011, 2012 and 2016, respectively. Data collected in 2011 and 2012 are referred to as Run 1 data, while data collected in 2016 are referred to as Run 2 data.

This paper is organised as follows. A brief description of the LHCb detector as well as the reconstruction and simulation software is given in Sect. 2. The selection of candidates is described in Sect. 3, and the first measurement of the branching fraction is presented in Sect. 4. An overview of the DP analysis formalism is given in Sect. 5. Details of the implementation of the DP fit are presented in Sect. 6. The evaluation of systematic uncertainties is given in Sect. 7. The results are summarised in Sect. 8.

Detector and simulation

The LHCb detector [21, 22] is a single-arm forward spectrometer covering the pseudorapidity range , designed for the study of particles containing or quarks. The detector includes a high-precision tracking system consisting of a silicon-strip vertex detector surrounding the pp interaction region, a large-area silicon-strip detector located upstream of a dipole magnet with a bending power of about , and three stations of silicon-strip detectors and straw drift tubes placed downstream of the magnet. The tracking system provides a measurement of the momentum, , of charged particles with a relative uncertainty that varies from 0.5% at low momentum to 1.0% at 200. The minimum distance of a track to a primary vertex (PV), the impact parameter, is measured with a resolution of , where is the component of the momentum transverse to the beam, in . Different types of charged hadrons are distinguished using information from two ring-imaging Cherenkov (RICH) detectors. Photons, electrons and hadrons are identified by a calorimeter system consisting of scintillating-pad and preshower detectors, an electromagnetic calorimeter and a hadronic calorimeter. Muons are identified by a system composed of alternating layers of iron and multiwire proportional chambers.

The online event selection is performed by a trigger [23], 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. At the hardware trigger stage, events are required to have a hadron with high transverse energy in the calorimeters. The software trigger requires a two-, three- or four-tracks secondary vertex with a significant displacement from any PV. At least one charged particle must have a large transverse momentum and be inconsistent with originating from a PV. A multivariate algorithm [24, 25] is used to identify secondary vertices that are consistent with b-hadron decays.

Simulated events, generated uniformly in the phase space of the or decay modes, are used to develop the selection, to validate the fit models and to evaluate the efficiencies entering the branching fraction measurement and the DP analysis. In the simulation, pp collisions are generated using Pythia  [26, 27] with a specific LHCb configuration [28]. Decays of hadronic particles are described by EvtGen  [29], in which final-state radiation is generated using Photos  [30]. The interaction of the generated particles with the detector, and its response, are implemented using the Geant4 toolkit [31, 32] as described in Ref. [33].

Selection

An initial offline selection comprising loose criteria is applied to reconstructed particles, where the associated trigger decision was due to the candidate. The final-state tracks are required to have , , and to be inconsistent with originating from any PV in the event. Loose particle identification (PID) criteria are applied, requiring the particles to be consistent with either the proton, kaon or pion hypothesis. All tracks are required to be within the acceptance of the RICH detectors (). Moreover, protons and antiprotons are required to have momenta larger than () to avoid kinematic regions where proton-kaon separation is limited for Run 1 (Run 2) data.

The candidates are required to have a small with respect to a PV, where is defined as the difference in the vertex-fit of a given PV reconstructed with and without the candidate under consideration. The PV providing the smallest value is associated to the candidate. The candidate is required to be consistent with originating from this PV by applying a criterion on the direction angle (DIRA) between the candidate momentum vector and the distance vector between the PV to the decay vertex. When building the candidates, the resolution on kinematic quantities such as the distribution, and the Dalitz variables that will be defined in Sect. 5, is improved by performing a kinematic fit [34] in which the candidate is constrained to originate from its associated PV, and its reconstructed invariant mass is constrained to the known mass [8].

A boosted decision tree (BDT) [35, 36] algorithm is used to further suppress the combinatorial background that arises when unrelated particles are combined to form a candidate. The training of the BDT is performed using simulated decays as the signal sample and candidates from the high-mass data sideband as the background sample, defined as the range. The input variables to the BDT classifiers are the same for Run 1 and 2 samples and comprise typical discriminating variables of b-hadron decays: the vertex-fit , , DIRA and flight distance significance of the reconstructed candidates; the maximum distance of closest approach between final-state particles; and the maximum and minimum p and of the proton and antiproton.

The requirements placed on the output of the BDT algorithm and PID variables are simultaneously optimised to maximise the figure of merit defined as . Here S is the observed yield before any BDT selection multiplied by the efficiency of the BDT requirement evaluated using simulated decays, while B is the the combinatorial background yield. The training of the BDT and the optimisation of the selection are performed separately for Run 1 and 2 data to accommodate for differences in the two data-taking periods.

Branching fraction measurement

The measurement of the branching fraction is performed relative to that of the normalisation channel, where the meson is also reconstructed in the decay mode. A two-stage fit procedure to the combined Run 1 and 2 data sample is used. In the first stage, an extended unbinned maximum-likelihood (UML) fit is performed to the distribution in order to separate the and background contributions. The RooFit package [37] is used to perform the fit, and the sPlot technique [38] is applied to assign weights for each candidate to subtract the background contributions. In the second stage, a weighted UML fit to the invariant-mass spectrum is performed to disentangle the , and nonresonant (NR) contributions. The efficiency-corrected yield ratio iswhere and are the observed and yields, while and are the total efficiencies, which are obtained from a combination of simulated and calibration samples. The branching fraction is determined aswhere , and are the external branching fractions taken from Ref. [8].

Signal and normalisation yields

The first-stage UML fit to the distribution is performed in the range. The signal decays, decays and various categories of background are present in this range. In addition to the combinatorial background, partially reconstructed backgrounds are present originating from b-hadron decays with additional particles that are not part of the reconstructed decay chain, such as a meson or a photon. Another source of background is b-hadron decays where one of the final-state particles has been incorrectly identified, which includes the decays and . The and decays are removed by excluding the mass range in the distribution and the range in the distribution, respectively. The latter veto also removes partially reconstructed b-hadron decays.

Both the and components are modelled by Hypatia functions [39]. The Hypatia distribution is a generalisation of the Crystall Ball function [40], where the Gaussian core of the latter is replaced by a hyperbolic core to take into account the distortion on the measured mass due to different sources of uncertainty. The Hypatia functions share a common resolution parameter, while the tail parameters are fixed to the values obtained from the corresponding simulated sample. The distributions of the misidentified and backgrounds are described by Crystal Ball functions, with parameters fixed to the values obtained from simulation. The combinatorial background is modelled using an exponential function. The masses of the and mesons, the resolution parameter of the Hypatia functions, the slope of the exponential function, and all the yields, are free to vary in the fit to the data. Using the information from the fit to the distribution, shown in Fig. 2, signal weights are computed and the background components are subtracted using the sPlot technique [38]. About decays are observed. Correlations between the and invariant-mass variables for both signal and background are found to be negligible.

Fig. 2

Distribution of the invariant mass. The solid blue curve is the projection of the total fit result. The components are shown in the legend

The second-stage UML fit is then performed to the weighted invariant-mass distribution in the mass range , which includes , and NR contributions. The invariant-mass distribution of candidates is described by the convolution of a nonrelativistic Breit–Wigner function and a Gaussian function describing resolution effects. Using simulated samples, the invariant-mass resolution is found to be . Given the width  [8], the impact of the detector resolution on the lineshape is small. The resonance, having a small natural width, is parametrised using an Hypatia function, with tail parameters fixed to the values obtained from the corresponding simulated sample. The same resolution parameter is used for the and contributions, which is free to vary in the fit to the data. The and masses are also floating, while the natural width is Gaussian constrained to the known value [8]. The NR contribution is parametrised with an exponential function with the slope free to vary in the fit. All yields are left unconstrained in the fit. A possible term describing the interference between the resonance and the NR S-wave is investigated and found to be negligible. The result of the fit to the weighted invariant-mass distribution is shown in Fig. 3. The yields of the and fit components, entering Eq. (1), are and , respectively.

Fig. 3

Distribution of the invariant mass in (left) linear and (right) logarithmic vertical-axis scale for weighted candidates obtained by using the sPlot technique. The solid blue curve is the projection of the total fit result. The full azure, tight-cross-hatched red and dashed-black line areas show the and NR contributions, respectively

Ratio of efficiencies

The ratio of efficiencies of Eq. (1) is obtained from and simulated samples, both selected using the same criteria used in data. Since these decays have the same final-state particles and similar kinematic distributions, the ratio of efficiencies is expected to be close to unity. The efficiencies are computed as the product of the geometrical acceptance of the LHCb detector, the reconstruction efficiency and the efficiency of the offline selection criteria, including the trigger and PID requirements. The efficiency of the PID requirements is obtained using calibration samples of pions, kaons and protons, as a function of the particle momentum, pseudorapidity and the multiplicity of the event, e.g. the number of charged particles in the event [41]. The final ratio of efficiencies is given bywhich is compatible with unity as expected.

Systematic uncertainties

Table 1 summarises the systematic uncertainties on the measurement of the ratio R of Eq. (1). Since the kinematic distributions of the signal and normalisation channel are similar, the uncertainties corresponding to the reconstruction and selection efficiencies largely cancel in the ratio of branching fractions. A new value of the ratio R is computed for each source of systematic uncertainty, and its difference with the nominal value is taken as the associated systematic uncertainty. The overall systematic uncertainty is assigned by combining all contributions in quadrature.

Table 1

Relative systematic uncertainties on the ratio R of Eq. (1). The total systematic uncertainty is obtained from the quadratic sum of the individual sources

Source	Systematic uncertainty (%)
Fixed shape parameters	0.8
Resolution model	0.3
NR pp¯ model	1.7
Efficiency ratio	1.1
Total	2.2

The systematic uncertainty arising from fixing the shape parameters of the Hypatia functions used to parametrise the and components is evaluated by repeating the fits and varying all shape parameters simultaneously. These shape parameters are varied according to normal distributions, taking into account the correlations between the parameters and with variances related to the size of the simulated samples.

To assign a systematic uncertainty arising from the model used to describe the detector resolution, the fits are repeated for each step replacing the Hypatia functions by Crystal Ball functions, whose parameters are obtained from simulation.

The systematic uncertainty associated to the parametrisation of the NR contribution is determined by replacing the exponential function with a linear function.

The systematic uncertainty associated to the determination of the efficiency involves contributions arising from the weighting procedure of the calibration samples used to determine the PID efficiencies. The granularity of the binning in the weighting procedure is halved and doubled.

The free shape parameters in the first stage UML fit lead to uncertainties that are not taken into account by the sPlot technique. In order to estimate this effect, these parameters are varied within their uncertainties and the signal weights are re-evaluated. The variations on the ratio R resulting from the second stage UML fit are found to be negligible.

Results

The ratio R is determined to bewhere the first uncertainty is statistical and the second systematic. The statistical uncertainty includes contributions from the per-candidate weights obtained using the sPlot technique. The value of R is used to compute the branching fraction using Eq. (2) which giveswhere the first uncertainty is statistical, the second systematic, and the third is due to the limited knowledge of the external branching fractions.

Dalitz plot formalism

The phase space for a three-body decay involving only pseudoscalar particles can be represented in a DP, where two of the three possible two-body invariant-mass-squared combinations, here and , are used to define the DP axes. However, given the sizeable natural width of the meson, the invariant mass is used instead of the known value of the mass [8] to compute the kinematic quantities such as and the helicity angles.

The isobar model [42-44] is used to write the decay amplitude as a coherent sum of amplitudes from resonant and NR intermediate processes aswhere are complex coefficients giving the relative contribution of each intermediate process. The complex functions describe the resonance dynamics and are normalised such that the integral of their squared magnitude over the DP is unityEach contribution is composed of the product of several factors. For a resonance, for instance,where is a normalisation constant and and are the momentum of the accompanying particle (the meson in this case) and the momentum of one of the resonance decay products, respectively, both evaluated in the rest frame. The X(z) terms are the Blatt–Weisskopf barrier factors [45] reported in Appendix A. The barrier radius, , is taken to be (corresponding to ) for all resonances. The term describes the angular probability distribution in the Zemach tensor formalism [46, 47], given by the equations reported in Appendix B. The function of Eq. (6) is the mass lineshape. Most of the resonant contributions are described by the relativistic Breit–Wigner (RBW) functionwhere the mass-dependent width is given byand is the value of for , being the pole mass of the resonance.

The amplitude parametrisations using RBW functions lead to unitarity violation within the isobar model if there are overlapping resonances or if there is a significant interference with a NR component, both in the same partial wave [48]. This is the case for the S-wave at low mass, where the resonance interferes strongly with a slowly varying NR S-wave component. Therefore, the S-wave at low mass is modelled using a modified LASS lineshape [49], given bywithand where and are the pole mass and width of the state, and a and r are the scattering length and the effective range, respectively. The parameters a and r depend on the production mechanism and hence on the decay under study. The slowly varying part (the first term in Eq. (9)) is not well modelled at high masses and it is set to zero for values above .

The probability density function for signal events across the DP, neglecting reconstruction effects, can be written aswhere the dependence of on the DP position has been suppressed for brevity. The natural width of the meson is set to zero when computing the DP normalisation shown in the denominator of Eq. (11). The effect of this simplification is determined when assessing the systematic uncertainties as described in Sect. 7.

The complex coefficients, given by in Eq. (4), depend on the choice of normalisation, phase convention and amplitude formalism. Fit fractions and interference fit fractions are convention-independent quantities that can be directly compared between different analyses. The fit fraction is defined as the integral of the amplitude for a single component squared divided by that of the coherent matrix element squared for the complete DP,In general, the fit fractions do not sum to unity due to the possible presence of net constructive or destructive interference over the whole DP area. This effect can be described by interference fit fractions defined for bywhere the dependence of and on the DP position is omitted.

Dalitz plot fit

The Laura package [50] is used to perform the unbinned DP fit, with the Run 1 and 2 subsamples fitted simultaneously using the JFIT framework [51]. The free parameters in the amplitude fit are in common between the two subsamples, while the signal and background yields and the maps describing the efficiency variations across the phase space, are different. Within the DP fit, the signal corresponds to decays, while the background comprises both combinatorial background and NR contributions. The likelihood function is given bywhere the index i runs over the candidates, k runs over the signal and background components, and is the yield of each component. The procedure to determine the signal and background yields is described in Sect. 6.1. The probability density function for the signal, , is given by Eq. (11) where the term is multiplied by the efficiency function described in Sect. 6.3. In order to avoid problems related to the imperfect parametrisation of the efficiencies at the DP borders, a veto of is applied around the DP, i.e. to the phase space boundaries of the , and distributions. This veto is used when determining the signal and background yields, and the probability density functions for the background, obtained as described in Sect. 6.2. The mass resolution is , which is much smaller than the meson width , the narrowest contribution to the DP; therefore, the resolution has negligible effects and is not considered further. The amplitude fits are repeated many times with randomised initial values to ensure the absolute minimum is found.

Signal and background yields

There is a non-negligible fraction of NR decays in the region of the meson. In order to separate the contributions of and NR decays, a two-dimensional (2D) UML fit to the and distributions is performed in the domain and . These ranges are chosen to avoid the misidentified decays reported in Sect. 4.1, and they also define the DP fit domain. The Run 1 and 2 2D mass fits are performed separately. The distributions of signal and NR decays are described by Hypatia functions. The distribution of the combinatorial background is parametrised using an exponential function. The distribution of signal decays is described by the same model described in Sect. 4.1. A possible component where genuine mesons are combined with random kaons and pions from the PV is investigated but found to be negligible. The meson mass, the resolution, the value of , the slopes of the exponential functions, and the yields, are free to vary in the 2D mass fits. The resolution and the meson natural width are Gaussian constrained to the value obtained in the fit to the weighted distribution of Sect. 4.1, and to the known value [8], respectively.

The yields of all fit components are reported in Table 2. Figure 4 shows the result of the 2D mass fits for the Run 1 and 2 subsamples that yield a total of approximately 2000 decays. The total yield of the component is lower than that reported in Sect. 4.1 since the fit ranges are reduced. The goodness of fit is validated using pseudoexperiments to determine the 2D pull, i.e. the difference between the fit model and data divided by the uncertainty.

Table 2

Yields of the components in the 2D mass fit to the joint [, ] distribution for the Run 1 and 2 subsamples

Component	Run 1	Run 2
	805±48	1065±56
(NR)	234±48	273±56
Combinatorial background	409±36	498±41

Fig. 4

Results of the 2D mass fit to the joint [, ] distribution for the a Run 1 projection, b Run 1 projection, c Run 2 projection, and d Run 2 projection. The legend is shown in the top left plot

Parametrisation of the backgrounds

The probability density functions for the combinatorial and NR background categories are obtained from the DP distribution of each background source, represented with a uniformly binned 2D histogram. In order to avoid artefacts related to the curved boundaries of the DP, the histograms are built in terms of the Square Dalitz plot (SDP) parametrised by the variables and which are defined in the range 0 to 1 and are given by where , are the kinematic boundaries of allowed in the decay, and is the helicity angle of the system (the angle between the and the mesons in the rest frame).

The combinatorial and NR background histograms are filled using the weights obtained by applying the sPlot technique to the joint [, ] distribution, merging the Run 1 and 2 data samples. Each histogram is scaled for the corresponding yield in the two subsamples. The combinatorial and NR background histograms for the Run 2 subsample are shown in Fig. 5. Statistical fluctuations in the histograms due to the limited size of the samples are smoothed by applying a 2D cubic spline interpolation.

Fig. 5

SDP distributions used in the DP fit to the Run 2 subsample for a combinatorial background and b NR background

The 2D mass fit described in Sect. 6.1 is repeated to the combined Run 1 and 2 data sample, and the sPlot technique is applied to determine the background-subtracted DP and SDP distributions shown in Fig. 6.

Fig. 6

Background-subtracted (top) DP and (bottom) SDP distributions corresponding to the total data sample used in the analysis. The structure corresponding to the resonance is evident. The veto of decays in the region is visible in the DP

Signal efficiency

Efficiency variation across the SDP is caused by the detector acceptance and by the trigger and offline selection requirements. The efficiency variation is evaluated with simulated samples generated uniformly across the SDP. Corrections are applied for known differences between data and simulation in PID efficiencies. The effect of the vetoes in the phase space is separately accounted for by the Laura package, setting to zero the signal efficiency within the vetoed regions. Therefore, the vetoes corresponding to the meson and the phase-space border are not applied when constructing the numerator of the efficiency histogram. The efficiency is studied separately for the Run 1 and 2 subsamples, and the resulting efficiency maps are shown in Fig. 7. Lower efficiency in regions with a low-momentum track is due to geometrical effects. Statistical fluctuations in the histograms due to the limited size of the simulated samples are smoothed by applying a 2D cubic spline interpolation.

Fig. 7

signal efficiency across the SDP for the a Run 1 and b Run 2 samples

Amplitude model with only contributions

In the absence of contributions from exotic resonances, only resonances are expected as intermediate states. The established mesons reported in Ref. [8] with , i.e. with masses within or slightly above the phase space boundary in decays, are used as a guide when building the model. Only those amplitudes providing significant improvements in the description of the data are included. This model is referred to as the baseline model and comprises the resonances shown in Table 3.

Table 3

Resonances included in the baseline model, where parameters and uncertainties are taken from Ref. [52]. The LASS lineshape also parametrise the S-wave in NR decays

Resonance	Mass [[]MeV]	Width [[]MeV]	J^P	Model
K^∗(892)^0	895.55±0.20	47.3±0.5	1^-	RBW
K^∗(1410)^0	1414±15	232±21	1^-	RBW
K0∗(1430)^0	1425±50	270±80	0^+	LASS
K2∗(1430)^0	1432.4±1.3	109±5	2^+	RBW
K^∗(1680)^0	1717±27	322±110	1^-	RBW
K0∗(1950)^0	1945±22	201±90	0^+	RBW

The S-wave at low mass is modelled with the LASS probability density function. The real and imaginary parts of the complex coefficients introduced in Eq. (4) are free parameters of the fit, except for the component, which is taken as the reference amplitude. Other free parameters in the fit are the scattering length (a) and the effective range (r) parameters of the LASS function, defined in Eq. (9). The mass and width of the meson are Gaussian constrained to the known values [8].

While it is possible to describe the and distributions well with contributions alone, the fit projection onto the distribution does not provide a good description of data, as shown in Fig. 8. In particular, a discrepancy around is evident.

Fig. 8

Projections of the data and amplitude fit using the baseline model onto a , c and e , with the same projections shown in b, d and f with a logarithmic vertical-axis scale. The veto of decays is visible in plot b. The S-wave component comprises the LASS and meson contributions. The components are described in the legend at the bottom

A variable is computed as a quantitative determination of the fit quality, using an adaptive 2D binning schema to obtain 144 equally populated bins. The baseline model yields a value, where ndof is the number of degrees of freedom. Including additional resonant states does not lead to significant improvements in the description of the data. These include established states such as the and mesons, the high mass resonance which falls outside the phase space limits, and the state which has not been seen in the final state thus far. The unestablished P-, D- and F-wave states predicted by the Godfrey–Isgur model [53] to decay into the final state were also tested.

Amplitude model with and contributions

A better description of the data is obtained by adding an exotic component to the contributions of Table 3. The resulting signal model consists of eight amplitudes: seven resonances and one NR term. The amplitudes are modelled in the same way as in the baseline model. Alternative models for the S-wave are used to assign systematic uncertainties as discussed in Sect. 7. In addition to the free parameters used in the baseline model, the isobar coefficients, mass and width of the resonance are left floating.

A likelihood-ratio test is used to discriminate between any pair of amplitude models based on the log-likelihood difference  [54]. Three quantum number hypotheses are probed for the resonance, repeating the amplitude fit for the assignments. The variations of the value with respect to the baseline model are , and 7.0, respectively. The model providing the best description of the data, referred to below as the nominal fit model, is obtained with the addition of a candidate with . The assignment is not considered further given the small variation in with respect to the additional four free parameters.

The LASS parameters obtained in the nominal fit model are , , and . The parameters of the candidate obtained in the nominal fit model are and . The values of the complex coefficients and fit fractions returned by the nominal fit model are shown in Table 4. The statistical uncertainties on all parameters of interest are calculated using large samples of simulated pseudoexperiments generated from the fit results in order to take into account the correlations between parameters and to guarantee the correct coverage of the uncertainties.

Table 4

Complex coefficients and fit fractions determined from the DP fit using the nominal model. Uncertainties are statistical only

Amplitude	Real part	Imaginary part	Fit fraction (%)
B^0→η_cK^∗(892)^0	1 (fixed)	0 (fixed)	51.4±1.9
B^0→η_cK^∗(1410)^0	0.17±0.07	0.11±0.08	2.1±1.1
B^0→η_cK^+π^- (NR)	-0.45±0.08	0.01±0.09	10.3±1.4
B^0→η_cK0∗(1430)^0	-0.62±0.09	-0.33±0.25	25.3±3.5
B^0→η_cK2∗(1430)^0	0.16±0.06	-0.23±0.05	4.1±1.5
B^0→η_cK^∗(1680)^0	-0.11±0.08	-0.18±0.06	2.2±2.0
B^0→η_cK0∗(1950)^0	0.27±0.04	0.04±0.14	3.8±1.8
B^0→Z_c(4100)^-K^+	-0.25±0.04	-0.01±0.08	3.3±1.1

Figure 9 shows the projections of the nominal fit model and the data onto , and invariant masses. A good agreement between the nominal fit model and the data is obtained. The value of the is for the nominal fit model. The fit quality is further discussed in Appendix C, where a comparison is reported of the unnormalised Legendre moments between data, the baseline and nominal models. The 2D pull distributions for the baseline and nominal models are reported as well.

Fig. 9

Projections of the data and amplitude fit using the nominal model onto a , c and e , with the same projections shown in b, d and f with a logarithmic vertical-axis scale. The veto of decays is visible in plot b. The S-wave component comprises the LASS and meson contributions. The components are described in the legend at the bottom

The significance of the candidate, referred to as the state in the following, is evaluated from the change in the likelihood of the fits with and without the component, assuming that this quantity, , follows a distribution with a number of degrees of freedom equal to twice the number of free parameters in its parametrisation [17, 55–57]. This assumption takes into account the look-elsewhere effect due to the floating mass and width of the . The validity of this assumption is verified using pseudoexperiments to predict the distribution of under the no- hypothesis, which is found to be well described by a probability density function with ndof 8. The statistical significance of the is in the nominal fit model. The quoted significance does not include the contribution from systematic uncertainties.

To discriminate between various assignments, fits are performed under alternative hypotheses. A lower limit on the significance of rejection of the hypothesis is determined from the change in the log-likelihood from the preferred hypothesis, assuming a distribution with one degree of freedom. The validity of this assumption is verified using pseudoexperiments to predict the distribution of under the disfavoured hypothesis. The statistical rejection of the hypothesis with respect to the hypothesis is .

Systematic effects must be taken into account to report the significance of the contribution and the discrimination of its quantum numbers. The fit variations producing the largest changes in the values of the mass, width or isobar coefficients of the exotic candidate are used to probe the sensitivity of the significance of the state to systematic effects, and to determine its quantum numbers, as described in Sect. 7.

Systematic uncertainties

Systematic uncertainties can be divided into two categories: experimental and model uncertainties. Among the experimental uncertainties, the largest changes in the values of the parameters of the candidate are due to the signal and background yields used in the amplitude fit, the SDP distributions of the background components, and the phase-space border veto applied on the parametrisation of the efficiencies. Among the model uncertainties, the largest effects are due to the treatment of the natural width of the meson within the DP fit and to the S-wave parametrisation. The DP fits using the baseline and nominal varied models are used to recompute the significance.

The signal and background yields used in the amplitude fit are fixed to the values obtained from the 2D mass fit. The statistical uncertainties on the yields are introduced into the amplitude fit by Gaussian constraining the yields within their statistical uncertainties and by repeating the fit.

The systematic uncertainties associated to the parametrisation of the background distributions are evaluated by varying the value in each bin within the statistical uncertainty prior to the spline interpolation. About 300 new background histograms are produced for both the combinatorial and NR background components. The resulting distribution follows a Gaussian distribution. The most pessimistic background parametrisation, corresponding to a value that is below of the Gaussian distribution, is considered when quoting the effect of this source on the significance of the state.

The phase-space border veto applied on the parametrisation of the efficiencies is removed to check the veto does not significantly affect the result.

The natural width of the meson is set to zero when computing the DP normalisations, calculated using the meson mass values resulting from the 2D UML fits described in Sect. 6.1. In order to associate a systematic uncertainty to the sizeable natural width, the amplitude fits are repeated computing the DP normalisations by using the and values, where and are the mass and natural width of the meson, respectively, obtained from the 2D UML fits.

The LASS model used to parametrise the low mass S-wave in the nominal fit is replaced with and resonances parametrised with RBW functions, and a NR S-wave component parametrised with a uniform amplitude within the DP.

The effect of the separate systematic sources to the significance of the are reported in Table 5. When including the most important systematic effect, corresponding to the pessimistic background parametrisation, the lowest significance for the candidate is given by . In order to evaluate the effect of possible correlated or anti-correlated sources of systematic uncertainty, the fits are repeated using the pessimistic background parametrisation together with the alternative S-wave model, and with mass values of the meson varied within the corresponding statistical uncertainty resulting from the 2D UML fit. The lower limit on the significance of the state is found to be .

Table 5

Significance of the contribution for the systematic effects producing the largest variations in the parameters of the candidate. The values obtained in the nominal amplitude fit are shown in the first row

Source	Δ(-2lnL)	Significance
Nominal fit	41.4	4.8σ
Fixed yields	45.8	5.2σ
Phase-space border veto	44.6	5.1σ
width	36.6	4.3σ
K^+π^- S-wave	31.8	3.9σ
Background	27.4	3.4σ

The discrimination between the and assignments is not significant when systematic uncertainties are taken into account, as reported in Table 6. When the S-wave model is varied, the two spin-parity hypotheses only differ by .

Table 6

Rejection level of the hypothesis with respect to the hypothesis for the systematic variations producing the largest variations in the parameters of the candidate. The values obtained in the nominal amplitude fit are shown in the first row

Source	Δ(-2lnL)	Significance
Default	18.6	4.3σ
Fixed yields	23.8	4.9σ
Phase-space border veto	24.4	4.9σ
width	4.2	2.0σ
Background	3.4	1.8σ
K^+π^- S-wave	1.4	1.2σ

Additional sources of systematic uncertainties are considered when evaluating the uncertainty on the mass and width of the resonance, and on the fit fractions obtained with the nominal model. These additional sources are: the efficiency variation across the SDP and a possible bias due to the fitting procedure, contributing to the experimental systematic uncertainties category; and the fixed parameters of the resonances in the amplitude model and the addition or removal of marginal amplitudes, contributing to the model systematic uncertainties category. For each source, the systematic uncertainty assigned to each quantity is taken as the difference between the value returned by the modified amplitude fit and nominal model fit result. The uncertainties due to all these sources are obtained by combining positive and negative deviations in quadrature separately.

The bin contents of the histograms describing the efficiency variation across the SDP are varied within their uncertainties prior to the spline interpolation, as is done for the systematic uncertainty associated to the background parametrisations. A possible source of systematic effects in the efficiency histograms is due to neighbouring bins varying in a correlated way. In order to evaluate this systematic uncertainty, 10 bins of the efficiency histograms are varied within their statistical uncertainty, and the neighbouring bins are varied by linear interpolation. The binning scheme of the control sample used to evaluate the PID performance is varied.

Pseudoexperiments are generated from the fit results using the nominal model in order to assign a systematic uncertainty due to possible amplitude fit bias.

Systematic uncertainties due to fixed parameters in the fit model are determined by repeating the fit and varying these parameters. The fixed masses and widths of the contributions are varied 100 times assigning a random number within the range defined by the corresponding uncertainties reported in Table 3. The Blatt–Weisskopf barrier radii, , are varied independently for and resonances between 3 and .

Systematic uncertainties are assigned from the changes in the results when the amplitudes due to the established and resonances, not contributing significantly in the baseline and nominal models, are included.

The total systematic uncertainties for the fit fractions are given together with the results in Sect. 8. The dominant experimental systematic uncertainty is due to either the phase-space border veto, related to the efficiency parametrisation, or the background distributions across the SDP, while the model uncertainties are dominated by the description of the S-wave.

The stability of the fit results is confirmed by several cross-checks. The addition of further high-mass states to the nominal model does not improve the quality of the fit. An additional amplitude decaying to is not significant, nor is an additional exotic amplitude decaying to . The meson resonant phase motion due to the sizeable natural width could affect the overall amplitude of Eq. (4), introducing interference effects with the NR contribution. In order to investigate this effect, the data sample is divided in two parts, containing candidates with masses below and above the meson peak, respectively. The results are compatible with those reported in Sect. 6 using the full data sample, supporting the argument that the effects due to the variation of the phase are negligible.

Table 7

Fit fractions and their uncertainties. The quoted uncertainties are statistical and systematic, respectively

Amplitude	Fit fraction (%)
B^0→η_cK^∗(892)^0	51.4±1.9-4.8+1.7
B^0→η_cK^∗(1410)^0	2.1±1.1-1.1+1.1
B^0→η_cK^+π^- (NR)	10.3±1.4-1.2+1.0
B^0→η_cK0∗(1430)^0	25.3±3.5-2.8+3.5
B^0→η_cK2∗(1430)^0	4.1±1.5-1.6+1.0
B^0→η_cK^∗(1680)^0	2.2±2.0-1.7+1.5
B^0→η_cK0∗(1950)^0	3.8±1.8-2.5+1.4
B^0→Z_c(4100)^-K^+	3.3±1.1-1.1+1.2

Results and summary

In summary, the first measurement of the branching fraction is reported and giveswhere the first uncertainty is statistical, the second systematic, and the third is due to limited knowledge of external branching fractions. The first Dalitz plot analysis of the decay is performed. A good description of data is obtained when including a charged charmonium-like resonance decaying to the final state with and . The fit fractions are reported in Table 7. The fit fractions for resonant and nonresonant contributions are converted into quasi-two-body branching fractions by multiplying by the branching fraction. The corresponding results are shown in Table 8. The branching fraction is compatible with the world-average value [8], taking into account the branching fraction. The values of the interference fit fractions are given in Table 9.

Table 8

Branching fraction results. The four quoted uncertainties are statistical, branching fraction systematic (not including the contribution from the uncertainty associated to the efficiency ratio, to avoid double counting the systematic uncertainty associated to the evaluation of the efficiencies), fit fraction systematic and external branching fractions uncertainties, respectively

Decay mode	Branching fraction (10^-5)
B^0→η_cK^∗(892)^0(→K^+π^-)	29.5±1.6±0.6-2.8+1.0±3.4
B^0→η_cK^∗(1410)^0(→K^+π^-)	1.20±0.63±0.02±0.63±0.14
B^0→η_cK^+π^- (NR)	5.90±0.84±0.11-0.69+0.57±0.68
B^0→η_cK0∗(1430)^0(→K^+π^-)	14.50±2.10±0.28-1.60+2.01±1.67
B^0→η_cK2∗(1430)^0(→K^+π^-)	2.35±0.87±0.05-0.92+0.57±0.27
B^0→η_cK^∗(1680)^0(→K^+π^-)	1.26±1.15±0.02-0.97+0.86±0.15
B^0→η_cK0∗(1950)^0(→K^+π^-)	2.18±1.04±0.04-1.43+0.80±0.25
B^0→Z_c(4100)^-K^+	1.89±0.64±0.04-0.63+0.69±0.22

Table 9

Symmetric matrix of the fit fractions (%) from the amplitude fit using the nominal model. The quoted uncertainties are statistical and systematic, respectively. The diagonal elements correspond to the values reported in Table 7

	K^∗(892)^0	K^∗(1410)^0	LASS NR	K0∗(1430)^0	K2∗(1430)^0	K^∗(1680)^0	K0∗(1950)^0	Z_c(4100)^-
K^∗(892)^0	51.4±1.9-4.8+1.7	1.7±1.9-1.4+2.4	0	0	0	-2.1±1.1-1.5+1.4	0	1.4±1.0-1.1+1.2
K^∗(1410)^0		2.1±1.1-1.1+1.1	0	0	0	-2.5±1.6-1.7+1.9	0	-0.4±0.4-0.5+0.7
LASS NR			10.3±1.4-1.2+1.0	-5.8±1.3-2.0+2.2	0	0	-3.2±2.8-1.4+4.9	1.11±0.23-0.35+0.54
K0∗(1430)^0				25.3±3.5-2.8+3.5	0	0	4.7±0.7-1.5+1.3	2.8±0.4-0.4+0.6
K2∗(1430)^0					4.1±1.5-1.6+1.0	0	0	0.00±0.31-0.26+0.76
K^∗(1680)^0						2.2±2.0-1.7+1.5	0	0.7±0.5-0.9+0.5
K0∗(1950)^0							3.8±1.8-2.5+1.4	0.6±0.5-1.1+0.8
Z_c(4100)^-								3.3±1.1-1.1+1.2

The significance of the candidate is more than three standard deviations when including systematic uncertainties. This is the first evidence for an exotic state decaying into two pseudoscalars. The favoured spin-parity assignments, and , cannot be discriminated once systematic uncertainties are taken into account, which prohibits unambiguously assigning the as one of the states foreseen by the models described in Sect. 1. Furthermore, the mass value of the charmonium-like state is above the open-charm threshold, in contrast with the predictions of such models. More data will be required to conclusively determine the nature of the candidate.

Comparison of the first four Legendre moments determined from background-subtracted data (black points) and from the results of the amplitude fit using the baseline model (red triangles) and nominal model (blue triangles) as a function of