Evidence for an ηc(1S)π-\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\eta _c(1S) \pi ^-$$\end{document} resonance in B0→ηc(1S)K+π-\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$B^0 \rightarrow \eta _c(1S) K^+\pi ^-$$\end{document} decays

A Dalitz plot analysis of B0→ηc(1S)K+π-\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${{B} ^0} \!\rightarrow \eta _c(1S) {{K} ^+} {{\pi } ^-} $$\end{document} decays is performed using data samples of pp collisions collected with the LHCb\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\text{ LHCb } $$\end{document} detector at centre-of-mass energies of s=7,8\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\sqrt{s}} =7,~8$$\end{document} and 13TeV\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$13{\,\mathrm {Te}\mathrm {V}} $$\end{document}, corresponding to a total integrated luminosity of 4.7fb-1\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$4.7 \,\text{ fb }^{-1} $$\end{document}. A satisfactory description of the data is obtained when including a contribution representing an exotic ηc(1S)π-\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\eta _c(1S) \pi ^-$$\end{document} resonant state. The significance of this exotic resonance is more than three standard deviations, while its mass and width are 4096±20-22+18MeV\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$4096 \pm 20~^{+18}_{-22} \,\mathrm {Me}\mathrm {V} $$\end{document} and 152±58-35+60MeV\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$152 \pm 58~^{+60}_{-35} \,\mathrm {Me}\mathrm {V} $$\end{document}, respectively. The spin-parity assignments JP=0+\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$J^P=0^+$$\end{document} and JP=1-\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$J^{P}=1^-$$\end{document} are both consistent with the data. In addition, the first measurement of the B0→ηc(1S)K+π-\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${{B} ^0} \!\rightarrow \eta _c(1S) {{K} ^+} {{\pi } ^-} $$\end{document} branching fraction is performed and gives B(B0→ηc(1S)K+π-)=(5.73±0.24±0.13±0.66)×10-4,\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \displaystyle \mathcal {B}({{B} ^0} \!\rightarrow \eta _c(1S) {{K} ^+} {{\pi } ^-} ) = (5.73 \pm 0.24 \pm 0.13 \pm 0.66) \times 10^{-4}, \end{aligned}$$\end{document}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][3][4][5][6][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][3][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 Z c (3900) − state, discovered by the BESIII collaboration in the J/ψ π − 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 Z c (3900) − state predicts an as-yet-unobserved charged charmoniumlike state with a mass of approximately [3800]MeV whose dominant decay mode is to the η c π − system. 2 Alternatively, states like the Z c (3900) − 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 η c π − system. The η c π − system carries isospin I = 1, G-parity G = −1, spin J = L and parity P = (−1) L , where L is the orbital angular momentum between the η c and the π − mesons. Lattice QCD calculations [14,15] predict the mass and quantum numbers of these states, comprising a I G (J P ) = 1 − (0 + ) state of mass [4025 ± 49]MeV , a I G (J P ) = 1 − (1 − ) state of mass [3770 ± 42]MeV , and a I G (J P ) = 1 − (2 + ) state of mass [4045 ± 44]MeV . The Z c (4430) − 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 η c π − 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 J P = 0 + candidate below the open-charm threshold that could decay into the η c π − final state. Therefore, the discovery of a charged charmonium-like meson in the η c π − sys- In this article, the B 0 → η c K + π − decay is studied for the first time, with the η c meson reconstructed using the p p decay mode. The decay is expected to proceed through K * 0 → K + π − intermediate states, where K * 0 refers to any neutral kaon resonance, following the diagram shown in Fig. 1a. If the decay also proceeds through exotic resonances in the η c π − system, denoted by Z − c states in the following, a diagram like that shown in Fig. 1b would contribute. The B 0 → η c K + π − 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.
The data sample used corresponds to an integrated luminosity of 4.7 fb −1 of pp collision data collected with the LHCb detector at centre-of-mass energies of √ s = 7, 8 and 13 TeV 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 B 0 → p pK + π − candidates is described in Sect. 3, and the first measurement of the B 0 → η c K + π − 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 2 < η < 5, designed for the study of particles containing b or c quarks. The detector includes a high-precision tracking system con-sisting 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 4 Tm, 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, p, of charged particles with a relative uncertainty that varies from 0.5% at low momentum to 1.0% at 200 GeV. The minimum distance of a track to a primary vertex (PV), the impact parameter, is measured with a resolution of (15 + 29/ p T ) µm, where p T is the component of the momentum transverse to the beam, in GeV. Different types of charged hadrons are distinguished using information from two ringimaging 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 B 0 → p pK + π − or B 0 → η c K + π − 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 B 0 candidate. The final-state tracks are required to have p > [1500]MeV , p T > [300]MeV , 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 (2.0 < η < 4.9). Moreover, protons and antiprotons are required to have momenta larger than [8]GeV ( [11]GeV ) to avoid kinematic regions where proton-kaon separation is limited for Run 1 (Run 2) data.
The B 0 candidates are required to have a small χ 2 IP with respect to a PV, where χ 2 IP is defined as the difference in the vertex-fit χ 2 of a given PV reconstructed with and without the candidate under consideration. The PV providing the smallest χ 2 IP value is associated to the B 0 candidate. The B 0 candidate is required to be consistent with originating from this PV by applying a criterion on the direction angle (DIRA) between the B 0 candidate momentum vector and the distance vector between the PV to the B 0 decay vertex. When building the B 0 candidates, the resolution on kinematic quantities such as the m( p p) distribution, and the Dalitz variables that will be defined in Sect. 5, is improved by performing a kinematic fit [34] in which the B 0 candidate is constrained to originate from its associated PV, and its reconstructed invariant mass is constrained to the known B 0 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 B 0 candidate. The training of the BDT is performed using simulated B 0 → p pK + π − decays as the signal sample and candidates from the high-mass data sideband as the background sample, defined as the 5450 < m( p pK + π − ) < 5550 MeV 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 χ 2 vtx , χ 2 IP , DIRA and flight distance significance of the reconstructed B 0 candidates; the maximum distance of closest approach between final-state particles; and the maximum and minimum p and p T 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 S/ √ S + B. Here S is the observed B 0 → p pK + π − 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 B 0 → η c K + π − branching fraction is performed relative to that of the B 0 → J/ψ K + π − normalisation channel, where the J/ψ meson is also reconstructed in the p p 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 m( p pK + π − ) distribution in order to separate the B 0 → p pK + π − 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 p p invariant-mass spectrum is performed to disentangle the η c , J/ψ , and nonresonant (NR) contributions. The efficiency-corrected yield ratio is where N η c and N J/ψ are the observed η c and J/ψ yields, while η c and J/ψ are the total efficiencies, which are obtained from a combination of simulated and calibration samples. The B 0 → η c K + π − branching fraction is determined as where

Signal and normalisation yields
The first-stage UML fit to the m( p pK + π − ) distribution is performed in the 5180 − 5430 MeV range. The B 0 → p pK + π − signal decays, B 0 s → p pK + π − 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  Fig. 2 Distribution of the p pK + π − invariant mass. The solid blue curve is the projection of the total fit result. The components are shown in the legend decays with additional particles that are not part of the reconstructed decay chain, such as a π 0 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 B 0 → p pπ + π − and B 0 s → p pK + K − . The D 0 → K + π − and − c → pK + π − decays are removed by excluding the mass range 1845 − 1885 MeV in the m(K + π − ) distribution and the range 2236−2336 MeV in the m( pK + π − ) distribution, respectively. The latter veto also removes partially reconstructed b-hadron decays.
Both the B 0 → p pK + π − and B 0 s → p pK + π − 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 B 0 → p pπ + π − and B 0 s → p pK + K − 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 B 0 and B 0 s 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 m( p pK + π − ) distribution, shown in Fig. 2, B 0 → p pK + π − signal weights are computed and the background components are subtracted using the sPlot technique [38]. About 3.0 × 10 4 B 0 decays are observed. Correlations between the p p and p pK + π − invariant-mass variables for both signal and background are found to be negligible.
The second-stage UML fit is then performed to the weighted p p invariant-mass distribution in the mass range 2700 − 3300 MeV, which includes η c , J/ψ , and NR B 0 → p pK + π − contributions. The p p invariant-mass distribution of η c candidates is described by the convolution of a nonrelativistic Breit-Wigner function and a Gaussian function describing resolution effects. Using simulated samples, the p p invariant-mass resolution is found to be ≈ [5]MeV . Given the width η c = [32.0 ± 0.8]MeV [8], the impact of the detector resolution on the η c lineshape is small. The J/ψ 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 η c and J/ψ contributions, which is free to vary in the fit to the data. The η c and J/ψ masses are also floating, while the η c natural width is Gaussian constrained to the known value [8]. The NR B 0 → p pK + π − 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 η c resonance and the NR p p S-wave is investigated and found to be negligible. The result of the fit to the weighted p p invariant-mass distribution is shown in Fig. 3. The yields of the B 0 → η c K + π − and B 0 → J/ψ K + π − fit components, entering Eq. (1), are 2105 ± 75 and 5899 ± 86, respectively.

Ratio of efficiencies
The ratio of efficiencies of Eq. (1) is obtained from B 0 → η c K + π − and B 0 → J/ψ K + π − 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 by which is compatible with unity as expected. 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 Candidates / (6 MeV)  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. The systematic uncertainty arising from fixing the shape parameters of the Hypatia functions used to parametrise the B 0 and J/ψ 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.

Systematic uncertainties
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 B 0 → p pK + π − 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 be where 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 B 0 → η c K + π − branching fraction using Eq. (2) which gives where 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 m 2 (K + π − ) and m 2 (η c π − ), are used to define the DP axes. However, given the sizeable natural width of the η c meson, the invariant mass m( p p) is used instead of the known value of the η c mass [8] to compute the kinematic quantities such as m 2 (η c K + ), m 2 (η c π − ) and the helicity angles.
The isobar model [42][43][44] is used to write the decay amplitude as a coherent sum of amplitudes from resonant and NR intermediate processes as where c j are complex coefficients giving the relative contribution of each intermediate process.
complex functions describe the resonance dynamics and are normalised such that the integral of their squared magnitude over the DP is unity where N is a normalisation constant and p and q are the momentum of the accompanying particle (the η c meson in this case) and the momentum of one of the resonance decay products, respectively, both evaluated in the K + π − rest frame. The X (z) terms are the Blatt-Weisskopf barrier factors [45] reported in Appendix A. The barrier radius, r BW , is taken to be [4]GeV −1 (corresponding to ∼ 0.8 fm) for all resonances. The Z ( p, q) term describes the angular probability distribution in the Zemach tensor formalism [46,47], given by the equations reported in Appendix B. The function T [m(K + π − )] of Eq. (6) is the mass lineshape. Most of the resonant contributions are described by the relativistic Breit-Wigner (RBW) function where the mass-dependent width is given by and q 0 is the value of | q| for m = m 0 , m 0 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 K + π − S-wave at low K + π − mass, where the K * 0 (1430) 0 resonance interferes strongly with a slowly varying NR S-wave component. Therefore, the K + π − S-wave at low mass is modelled using a modified LASS lineshape [49], given by with cot δ B = 1 a| q| and where m 0 and 0 are the pole mass and width of the K * 0 (1430) 0 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 The probability density function for signal events across the DP, neglecting reconstruction effects, can be written as where the dependence of A on the DP position has been suppressed for brevity. The natural width of the η c 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 c j 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 i < j by where the dependence of F ( * ) i and A 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 B 0 → η c K + π − decays, while the background comprises both combinatorial background and NR B 0 → p pK + π − contributions. The likelihood function is given by where the index i runs over the N c candidates, k runs over the signal and background components, and N k 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, P sig , is given by Eq. (11) where the |A[m 2 (K + π − ), m 2 (η c π − )]| 2 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 ± [70]MeV is applied around the DP, i.e. to the phase space boundaries of the m(K + π − ), m(η c π − ) and m(η c K + ) 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 K + π − mass resolution is ≈ [5]MeV , which is much smaller than the K * (892) 0 meson width K * (892) 0 ≈ [50]MeV , 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 B 0 → p pK + π − decays in the region of the η c meson. In order to separate the contributions of B 0 → η c K + π − and NR B 0 → p pK + π − decays, a two-dimensional (2D) UML fit to the m( p pK + π − ) and m( p p) distributions is performed in the domain 5220 < m( p pK + π − ) < 5340 MeV and 2908 < m( p p) < 3058 MeV. 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 m( p pK + π − ) distributions of B 0 → η c K + π − signal and NR B 0 → p pK + π − decays are described by Hypatia functions. The m( p pK + π − ) distribution of the combinatorial background is parametrised using an exponential function. The m( p p) distribution of B 0 → η c K + π − signal decays is described by the same model described in Sect  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 B 0 → η c K + π − decays. The total yield of the B 0 → η c K + π − 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.

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 m and θ which are defined in the range 0 to 1 and are given by where m max K + π − = m B 0 − m η c , m min K + π − = m K + + m π − are the kinematic boundaries of m(K + π − ) allowed in the B 0 → η c K + π − decay, and θ(K + π − ) is the helicity angle of the K + π − system (the angle between the K + and the η c mesons in the K + π − rest frame).
The combinatorial and NR background histograms are filled using the weights obtained by applying the sPlot technique to the joint [m( p pK + π − ), m( p p)] 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.
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.

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 D 0 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.

Amplitude model with only K + π − contributions
In the absence of contributions from exotic resonances, only K + π − resonances are expected as intermediate states. The established K * 0 → K + π − mesons reported in Ref. [8] with m(K * 0 ) m(B 0 ) − m(η c ), i.e. with masses within or slightly above the phase space boundary in B 0 → η c K + π − 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. The S-wave at low K + π − mass is modelled with the LASS probability density function. The real and imaginary parts of the complex coefficients c j introduced in Eq. (4) are free parameters of the fit, except for the K * (892) 0 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 K * 0 (1430) 0 meson are Gaussian constrained to the known values [8].
While it is possible to describe the m(K + π − ) and m(η c K + ) distributions well with K + π − contributions alone, the fit projection onto the m(η c π − ) distribution does not provide a good description of data, as shown in Fig. 8. In particular, a discrepancy around m(η c π − ) ≈ [4.1]GeV is evident.
A χ 2 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 χ 2 /ndof = 195/129 = 1.5 value, where ndof is the number of degrees of freedom. Including additional K + π − resonant states does not lead to significant improvements in the description of the data. These include established states such as the K * 3 (1780) 0 and K * 4 (2045) 0 mesons, the high mass K * 5 (2380) 0 resonance which falls outside the phase space limits, and the K * 2 (1980) 0 state which has not been seen in the K + π − final state thus far. The unestablished P-, Dand F-wave K + π − states predicted by the Godfrey-Isgur model [53] to decay into the K + π − final state were also tested.
6.5 Amplitude model with K + π − and η c π − contributions A better description of the data is obtained by adding an exotic Z − c → η c π − component to the K + π − contributions of Table 3. The resulting signal model consists of eight amplitudes: seven resonances and one NR term. The K + π − amplitudes are modelled in the same way as in the baseline model. Alternative models for the K + π − 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 Z − c resonance are left floating.
A likelihood-ratio test is used to discriminate between any pair of amplitude models based on the log-likelihood difference (−2 ln L) [54]. Three quantum number hypotheses are probed for the Z − c resonance, repeating the amplitude fit for the J P = 0 + , 1 − and 2 + assignments.  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. Figure 9 shows the projections of the nominal fit model and the data onto m(K + π − ), m(η c π − ) and m(η c K + ) invariant masses. A good agreement between the nominal fit model and the data is obtained. The value of the χ 2 /ndof is 164/125 = 1.3 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.
The significance of the Z − c candidate, referred to as the Z c (4100) − state in the following, is evaluated from the change in the likelihood of the fits with and without the Z c (4100) − component, assuming that this quantity, (−2 ln L), follows a χ 2 distribution with a number of degrees of freedom equal to twice the number of free parameters in its parametrisation [17,[55][56][57]. This assumption takes into account the look-elsewhere effect due to the floating mass and width of the Z c (4100) − . The validity of this assumption is verified using pseudoexperiments to pre- The veto of B 0 → p p D 0 decays is visible in plot b. The K + π − S-wave component comprises the LASS and K * 0 (1950) 0 meson contributions. The components are described in the legend at the bottom hypothesis, which is found to be well described by a χ 2 probability density function with ndof = 8. The statistical significance of the Z c (4100) − is 4.8σ in the nominal fit model. The quoted significance does not include the contribution from systematic uncertainties.
To discriminate between various J P assignments, fits are performed under alternative J P hypotheses. A lower limit on the significance of rejection of the J P = 0 + hypothesis is determined from the change in the log-likelihood from the preferred hypothesis, assuming a χ 2 distribution with one degree of freedom. The validity of this assumption is verified using pseudoexperiments to predict the distribution of (−2 ln L) under the disfavoured J P = 0 + hypothesis. The statistical rejection of the J P = 0 + hypothesis with respect to the J P = 1 − hypothesis is 4.3σ .
Systematic effects must be taken into account to report the significance of the Z c (4100) − 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 Z c (4100) − 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 Z c (4100) − candidate are due to the signal and background yields used in the amplitude fit, the SDP distributions of the background components, and the phasespace 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 η c meson within the DP fit and to the K + π − 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 (−2 ln L) distribution follows a Gaussian distribution. The most pessimistic background parametrisation, corresponding to a (−2 ln L) value that is below 3σ of the Gaussian distribution, is considered when quoting the effect of this source on the significance of the Z c (4100) − 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 η c meson is set to zero when computing the DP normalisations, calculated using the η c meson mass values resulting from the 2D UML fits described in Sect. 6.1. In order to associate a systematic uncertainty to the sizeable η c natural width, the amplitude fits are repeated computing the DP normalisations by using the m η c + η c and m η c − η c values, where m η c and η c are the mass and natural width of the η c meson, respectively, obtained from the 2D UML fits. The LASS model used to parametrise the low mass K + π − S-wave in the nominal fit is replaced with K * 0 (1430) 0 and K * 0 (700) 0 resonances parametrised with RBW functions, and a NR S-wave K + π − component parametrised with a uniform amplitude within the DP.
The effect of the separate systematic sources to the significance of the Z c (4100) − are reported in Table 5. When including the most important systematic effect, corresponding to the pessimistic background parametrisation, the lowest significance for the Z c (4100) − candidate is given by 3.4σ . 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 K + π − S-wave model, and with mass values of the η c meson varied within the corresponding statistical uncertainty resulting from the 2D UML fit. The lower limit on the significance of the Z c (4100) − state is found to be 3.2σ .
The discrimination between the J P = 0 + and J P = 1 − assignments is not significant when systematic uncertainties are taken into account, as reported in Table 6. When the Swave model is varied, the two spin-parity hypotheses only differ by 1.2σ .
Additional sources of systematic uncertainties are considered when evaluating the uncertainty on the mass and width of the Z c (4100) − 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 K + π − 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, r BW , are varied independently for K + π − and η c π − resonances between 3 and [5]GeV −1 .
Systematic uncertainties are assigned from the changes in the results when the amplitudes due to the established K * 3 (1780) 0 and K * 4 (2045) 0 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 phasespace 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 K + π − S-wave.
The stability of the fit results is confirmed by several crosschecks. The addition of further high-mass K * 0 states to the nominal model does not improve the quality of the fit. An additional amplitude decaying to η c π − is not significant, nor is an additional exotic amplitude decaying to η c K + . The  Table 8 Branching fraction results. The four quoted uncertainties are statistical, B 0 → η c K + π − 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 ) η c meson resonant phase motion due to the sizeable natural width could affect the overall amplitude of Eq. (4), introducing interference effects with the NR pp contribution. In order to investigate this effect, the data sample is divided in two parts, containing candidates with masses below and above the η c 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 η c phase are negligible.

Results and summary
In summary, the first measurement of the B 0 → η c K + π − branching fraction is reported and gives where 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 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  Table 7. The fit fractions for resonant and nonresonant contributions are converted into quasi-two-body branching fractions by multiplying by the B 0 → η c K + π − branching fraction. The corresponding results are shown in Table 8. The B 0 → η c K * (892) 0 branching fraction is compatible with the world-average value [8], taking into account the K * (892) 0 → K + π − branching fraction. The values of the interference fit fractions are given in Table 9.
The significance of the Z c (4100) − 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, J P = 0 + and J P = 1 − , cannot be discriminated once systematic uncertainties are taken into account, which prohibits unambiguously assigning the Z c (4100) − as one of the states foreseen by the models described in Sect. 1. Furthermore, the mass value of the Z c (4100) − 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 Z c (4100) − candidate.

Appendix: A Blatt-Weisskopf barrier factors
The Blatt-Weisskopf barrier factors [45] X (z), where z = | q|r BW or | p|r BW with r BW being the barrier radius, are given by where z 0 is the value of z when the invariant mass is equal to the pole mass of the resonance and L is the orbital angular momentum between the resonance children. Since the latter are scalars, L is equal to the spin of the resonance. Since the B 0 meson and the accompanying particle in the decay are scalars as well, L is also equal to the orbital angular momentum between the resonance and the accompanying particle in the decay.

B: Angular probability distributions
Using the Zemach tensor formalism [46,47], the angular probability distributions Z ( p, q) are given by Comparisons of the first four Legendre moments determined from background-subtracted data and from the amplitude fit results using the baseline and nominal model are reported in Figs. 10, 11 and 12 for the m(K + π − ), m(η c π − ) and m(η c K + ) projections, respectively. The 2D pull distributions for the baseline and nominal models are reported in Figs. 13 and 14