Study of the lineshape of the $\chi_{c1}(3872)$ state

A study of the lineshape of the $\chi_{c1}(3872)$ state is made using a data sample corresponding to an integrated luminosity of $3\,$fb$^{-1}$ collected in $pp$ collisions at centre-of-mass energies of $7$ and $8\,$TeV with the LHCb detector. Candidate $\chi_{c1}(3872)$ mesons from $b$-hadron decays are selected in the $ J/\psi \pi^+ \pi^-$ decay mode. Describing the lineshape with a Breit--Wigner function, the mass splitting between the $\chi_{c1}(3872)$ and $\psi(2S)$ states, $\Delta m$, and the width of the $\chi_{c1}(3872)$ state, $\Gamma_{\mathrm{BW}}$, are determined to be \begin{eqnarray*} \Delta m&=&185.588 \pm 0.067 \pm 0.068\,{\mathrm{MeV}} \,, \\ \Gamma_{\mathrm{BW}}&=&\phantom{00}1.39\phantom{0} \pm 0.24\phantom{0} \pm 0.10\phantom{0}\,{\mathrm{MeV}} \,, \end{eqnarray*} where the first uncertainty is statistical and the second systematic. Using a Flatt\'e-inspired lineshape, two poles for the $~\chi_{c1}(3872)~$state in the complex energy plane are found. The dominant pole is compatible with a quasi-bound $D^0\bar{D}^{*0}~$state but a quasi-virtual state is still allowed at the level of $2~$standard deviations.

A striking feature of the χ c1 (3872) state is the proximity of its mass to the sum of the D * 0 and D 0 meson masses. Accounting for correlated uncertainties due to the knowledge of the kaon mass, this sum is evaluated to be m D 0 + m D * 0 = 3871.70 ± 0.11 MeV. 1 The molecular interpretation of the χ c1 (3872) state requires it to be a bound state. Assuming a Breit-Wigner lineshape, this implies that δE ≡ m D 0 + m D * 0 − m χ c1 (3872) > 0. Current knowledge of δE is limited by the uncertainty on the χ c1 (3872) mass, motivating a more precise determination of this quantity. The nature of the χ c1 (3872) state can also be elucidated by studies of its lineshape. This has been analysed by several experiments assuming a Breit-Wigner function [3,5,13]. The current upper limit on the natural width, Γ BW , is 1.2 MeV at 90% confidence level [14].
In this analysis a sample of χ c1 (3872) → J/ψπ + π − candidates produced in inclusive b-hadron decays is used to measure precisely the mass and to determine the lineshape of the χ c1 (3872) meson. Studies are made assuming both a Breit-Wigner lineshape and a Flatté-inspired model that accounts for the opening up of the D 0 D * 0 threshold [15,16]. The analysis uses a data sample corresponding to an integrated luminosity of 3 fb −1 of data collected in pp collisions at centre-of-mass energies of 7 TeV and 8 TeV during 2011 and 2012 using the LHCb detector.

Detector and simulation
The LHCb detector [17,18] 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 consisting of a silicon-strip vertex detector surrounding the pp interaction region [19], a large-area silicon-strip detector (TT) 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 [20] placed downstream of the magnet. The tracking system provides a measurement of momentum, p, of charged particles with a relative uncertainty that varies from 0.5% at low momentum to 1.0% at 200 GeV. As described in Refs. [21,22] the momentum scale is calibrated using samples of J/ψ → µ + µ − and B + → J/ψ K + decays collected concurrently with the data sample used for this analysis. The relative accuracy of this procedure is estimated to be 3 × 10 −4 using samples of other fully reconstructed b hadrons, Υ and K 0 S mesons. The minimum distance of a track to a primary vertex (PV), the impact parameter (IP), 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 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 [23].
The online event selection is performed by a trigger [24], which consists of a hardware stage, based on information from the calorimeter and muon systems, followed by a software stage, where a full event reconstruction is made. Candidate events are required to pass the hardware trigger, which selects muon and dimuon candidates with high p T based upon muon system information. The subsequent software trigger is composed of two stages. The first performs a partial event reconstruction and requires events to have two well-identified oppositely charged muons and that the mass of the pair is larger than 2.7 GeV. The second stage performs a full event reconstruction. Events are retained for further processing if they contain a displaced µ + µ − vertex. The decay vertex is required to be well separated from each reconstructed PV of the proton-proton interaction by requiring the distance between the PV and the µ + µ − vertex divided by its uncertainty to be greater than three.
To study the properties of the signal and the most important backgrounds, simulated sampels of pp collisions are generated using Pythia [25] with a specific LHCb configuration [26]. Decays of hadronic particles are described by EvtGen [27], in which final-state radiation is generated using Photos [28]. The interaction of the generated particles with the detector, and its response, are implemented using the Geant4 toolkit [29] as described in Ref. [30]. For the study of the lineshape it is important that the simulation models well the mass resolution. The simulation used in this study reproduces the observed mass resolution for selected samples of B + → J/ψ K + , B 0 → J/ψ K + π − , B 0 s → J/ψ φ and B + → J/ψ K + π + π − decays within 5%. To further improve the agreement for the mass resolution between data and simulation, scale factors are determined using a large sample of ψ(2S) → J/ψ π + π − decays collected concurrently with the χ c1 (3872) sample. This is discussed in detail in Sec. 5.

Selection
The selection of χ c1 (3872) → J/ψ π + π − candidates from b-hadron decays is performed in two steps. First, loose selection criteria are applied that reduce the background from random combinations of tracks significantly while retaining high signal efficiency. Subsequently, a multivariate selection is used to further reduce this combinatorial background. In both steps, the selection requirements are chosen to reduce background whilst selecting well reconstructed candidates. The requirements are optimised using simulated signal decays together with a sample of selected candidates in data where the charged pions have the same sign. The latter sample is found to be a good proxy to describe the background shape. Though the selection criteria are tuned using the χ c1 (3872) simulation sample, the ψ(2S) → J/ψ π + π − decay mode is also selected with high efficiency and used for calibration.
The selection starts from a pair of oppositely charged particles, identified as muons. Incorrectly reconstructed tracks are suppressed by imposing a requirement on the output of a neural network trained to discriminate between these and trajectories from real particles. To select J/ψ → µ + µ − candidates, the two muons are required to originate from a common vertex that is significantly displaced from any PV. The difference between the reconstructed invariant mass of the pair and the known value of the J/ψ mass [31] is required to be within three times the uncertainty on the reconstructed mass of the µ + µ − pair.
Pion candidates are selected using the same track-quality requirements as the muons. Information from the muon system is used to reject pions that decayed in the spectrometer since these pions tend to have poorly reconstructed trajectories which result in χ c1 (3872) candidates with worse mass resolution. Combinatorial background is suppressed by requiring that the χ 2 IP of the pion candidates, defined as the difference between the χ 2 of the PV reconstructed with and without the considered particle, is larger than four for all PVs. Good pion identification is ensured by applying a requirement on a variable that combines information from the RICH detectors with kinematic and track quality information. Since the pions produced in χ c1 (3872) decays have relatively small p T , only a loose requirement on the transverse momentum (p T > 200 MeV) is imposed. In addition, the pion candidates are required to have p < 50 GeV. This requirement rejects candidates with poor momentum resolution and has an efficiency of 99.5 %.
To create χ c1 (3872) candidates, J/ψ candidates are combined with pairs of oppositely charged pions. To improve the mass resolution a kinematic vertex fit [32] is made. constraining the J/ψ invariant mass to its known value [31]. The reduced χ 2 of the fit, χ 2 fit /ndf, is required to be less than five. Candidates with a mass uncertainty greater than 5.0 MeV are rejected. Finally, requiring the Q-value of the decay to be below 200 MeV substantially reduces the background whilst retaining 96% of the χ c1 (3872) signal. Here the Q-value is defined as and m π + π − are the reconstructed masses of the final state combinations.
The final step of the selection process is based on a neural network classifier [33][34][35][36]. This is trained on a simulated sample of inclusive b → χ c1 (3872)X decays and the samesign pion sample in data. Simulated samples are corrected to reproduce kinematical distributions of the ψ(2S) mesons observed on data. The training is performed separately for the 2011 and 2012 data samples. Twelve variables that give good separation between signal and background are considered: the pseudorapidity and transverse momentum of the two pion candidates, the χ 2 IP for each of the two pions, the pseudorapidity and transverse momentum of the χ c1 (3872) candidate, the χ 2 of the two-track vertex fit for the pions, the χ 2 fit /ndf, the flight distance χ 2 of the candidate calculated using the reconstructed primary and secondary vertices, and the total number of hits in the TT detector. All these variables show good agreement between simulation and data. The optimal cut on the classifier output is chosen using pseudoexperiments so as to minimise the uncertainty on the measured χ c1 (3872) mass.

Mass model
The observed invariant mass distribution of the J/ψ π + π − system, m J/ψ π + π − , for the ψ(2S) and χ c1 (3872) resonances is a convolution of the natural lineshape with the detector resolution. For the ψ(2S) resonance the lineshape is well described by a Breit-Wigner function. The situation for the χ c1 (3872) meson is more complex. Previous measurements have assumed a Breit-Wigner resonance shape. However, as discussed in Refs. [12,15,16], this is not well motivated due to the proximity of the D 0 D * 0 threshold. Several other alternative lineshapes have been proposed in the literature [15,16,37,38]. In this analysis two lineshapes for the χ c1 (3872) meson are considered in detail, a Breit-Wigner and a Flatté-inspired model [15,16]. These models are investigated in the next sections. The S-wave threshold resonance model described in Ref. [37,38], that accounts for the non-zero width of the D * 0 meson, was considered but did not fit the data well. If the mass is close to the D 0 D * 0 threshold, this model is not able to accommodate a value of the natural width much larger than Γ D * 0 = 65.5 ± 15.4 keV [37]. As will be discussed below, the study presented here favours larger values of the natural width.
The analysis proceeds in two steps. First, unbinned maximum-likelihood fits are made to the m J/ψ π + π − distribution in the region around the ψ(2S) mass (Sec. 5). These measured values of the ψ(2S) mass and mass resolution are used to control systematic uncertainties in the subsequent fits to the m J/ψ π + π − distribution in the χ c1 (3872) mass region. For both sets of fits the natural lineshape is convolved with a resolution model developed using the simulation. The application of the J/ψ mass constraint in the fit [32] results in the mass resolution being dominated by the kinematics of the pion pair. In particular, the resolution is worse for higher values of the total momentum of the pion pair, p π + π − . Consequently, the analysis is performed in three p π + π − bins chosen to contain an approximately equal number of signal candidates: p π + π − < 12 GeV, 12 ≤ p π + π − < 20 GeV and 20 ≤ p π + π − < 50 GeV. The core mass resolution for the χ c1 (3872) state varies monotonically between 2.4 MeV and 3.0 MeV between the lowest-p π + π − and highest-p π + π − bin. Possible differences in data taking conditions are allowed for by dividing the data according to the year of collection resulting in a total of six data samples.
The resolution model is studied using simulation. In each p π + π − bin the mass resolution is modelled with the sum of a narrow Crystal Ball function [39] combined with a wider Gaussian function. The Crystal Ball function has a Gaussian core and two parameters that describe the power-law tail. The simulation is also used to determine the value of the transition point between the core and the power law tail, a, as a multiple of the width, σ, of the Gaussian core. The value of the exponent of the power law, n, is allowed to vary in the data fits with a Gaussian constraint to the value obtained in the simulation applied. When fitting the χ c1 (3872) mass region in data the values of the core resolution, σ, for the Gaussian and Crystal Ball functions are taken from simulation up to an overall scale factor, s f , that accounts for residual discrepancies between data and simulation. For each p π + π − data sample the value of s f is determined in the corresponding fit to the ψ(2S) mass region and applied as a Gaussian constraint. The systematic uncertainty associated with the choice of the signal model is assessed by replacing the nominal model with the sum of either two Crystal Ball or Gaussian functions.
The shape of the combinatorial background is studied using the same-sign data sample as well as samples of simulated inclusive b → J/ψX decays. Based upon these studies, the background is modelled by the form (m J/ψ π + π − − m J/ψ − 2m π ± ) c 0 e −m J/ψ π + π − /c 1 , where c 0 is fixed to 3.6 based on fits to the same-sign data. Variations of this functional form together with other models (e.g. exponential or polynomial functions) are used as systematic variations. In total, seven different background forms are considered.

ψ(2S) mass
Since the ψ(2S) state is narrow and away from the phase-space limits, a spin-0 relativistic Breit-Wigner function is used to model the lineshape. A spin-1 Breit-Wigner function is considered as part of the systematic uncertainties and found to give identical results. This lineshape is convolved with the default resolution model described in Sec. 4 and a fit to J/ψ π + π − mass is performed in each of the six p π + π − data samples. The natural width of the ψ(2S) is fixed to the known value [31]. Figure 1 shows the m J/ψ π + π − distributions and fit projections for each data sample and Table 1 summarises the resulting parameters of interest. Binning the data and calculating the χ 2 probability of consistency with the fit model gives the values greater than 5% for all fits. The fitted values of the ψ(2S) mass agree with the known value [31] within the uncertainty of the calibration procedure. The values of s f are consistent with the expectation that the simulation reproduces the mass resolution in data at the level of 5% or better. When applied as constraint in the fit to the χ c1 (3872) region, additional uncertainties on s f are considered. Accounting for the finite size of the simulation samples, the background modelling and the assumption that the ψ(2S) calibration factor can be applied to the χ c1 (3872) candidates, the uncertainty on s f is 0.02, independent of the bin. The values of s f in Table 1 are applied as Gaussian constraints in the fits to the χ c1 (3872) region with an uncertainty of 0.02.

Breit-Wigner mass and width of the χ c1 (3872) state
To extract the Breit-Wigner lineshape parameters of the χ c1 (3872) meson, a fit is made to the mass range 3832 < m J/ψ π + π − < 3912 MeV in each of the six p π + π − data samples described in Sec. 4. A spin-0 relativistic Breit-Wigner is used, as in Ref. [9], For each data sample the mass difference between the ψ(2S) and χ c1 (3872) meson, ∆m, is measured relative to the measured mass of the ψ(2S) state rather than the absolute mass. This minimises the systematic uncertainty due to the momentum scale. The fit in each bin has seven free parameters: ∆m, the natural width Γ BW , the background parameter c 1 , the resolution scale factor s f , the tail parameter n, and the signal and background yields. As in Sec. 5, a Gaussian constraint is applied to n based on the simulation. The parameter s f is constrained to the result of Sec. 5. The fit procedure  Year   is validated using both the simulation and pseudoexperiments. No significant bias is found and the uncertainties estimated by the fit agree with the spread observed in the pseudoexperiments. These studies show that, values of Γ BW larger than 0.6 MeV can reliably be determined.
For the six p π + π − data samples the J/ψ π + π − mass distributions in the χ c1 (3872) region and fits are shown in Fig. 2 LHCb 2012 2011 χ c1 (3872) background total and calculating the χ 2 probability of consistency with the fit model gives values much larger than 5% for all bins apart from the high-momentum bin in the 2012 data where the probability is 2%. The values of ∆m and Γ BW are consistent between the bins giving confidence in the results.
A simultaneous fit is made to the six data samples with ∆m and Γ BW as shared parameters. This gives ∆m = 185.588 ± 0.067 MeV and Γ BW = 1.39 ± 0.24 MeV, where Year The dominant systematic uncertainty on the mass difference ∆m arises from the 3 × 10 −4 relative uncertainty on the momentum scale. Its effect is evaluated by adjusting the four-vectors of the pions by this amount and repeating the analysis. The bias on ∆m from QED radiative corrections is determined to be (−10 ± 14) keV using the simulation, which uses Photos [28] to model this effect. The measured value of ∆m is corrected by this value and the uncertainty considered as a systematic error. The small uncertainty on the fitted values of the ψ(2S) mass is also propagated to the ∆m value. Biases arising from the modelling of the resolution and the treatment of the background shape are evaluated to be 2 keV using the discrete profiling method described in Ref. [40]. The uncertainties on the ∆m measurement are summarised in Table 3. Combining all uncertainties, the mass splitting between the χ c1 (3832) and ψ(2S) mesons is determined as where the first uncertainty is statistical and the second is systematic. The value of ∆m can be translated into an absolute measurement of the χ c1 (3872) mass using m ψ(2S) = 3686.097 ± 0.010 MeV from Ref. [31], yielding where the third uncertainty is due to the knowledge of the ψ(2S) mass. For these measurements it is assumed that interference effects with other partially reconstructed b-hadron decays do not affect the lineshape. This assumption is reasonable since many exclusive b-hadron decays contribute to the final sample, and the χ c1 (3872) state is narrow. This assumption has been explored in pseudoexperiments varying the composition and phases of the possible decay amplitudes that are likely to contribute to the observed data set. These studies conservatively limit the size of any possible effect on m χ c1 (3872) to be less than 40 keV.
Systematic uncertainties on Γ BW largely arise from the modelling of the resolution function. The uncertainties on s f and n are propagated to Γ BW via Gaussian constraints and hence are included in the statistical uncertainty, to which they contribute 0.05 MeV. The uncertainty due to the choice of the signal and background models is evaluated using the discrete profiling method with the alternative models described in Sec. 4. Based upon these studies an uncertainty of 0.10 MeV is assigned. The uncertainty due to possible differences in the p T distribution between data and simulation is evaluated by weighting the simulation to achieve better agreement and lead to a 0.01 MeV uncertainty. Summing these values in quadrature gives a total uncertainty of 0.1 MeV. The value of Γ BW , including systematic uncertainties, differs from zero by more than 5 standard deviations. Fits were also made fixing Γ BW to zero and allowing s f to float in each bin without constraint. The value of s f obtained is between 1.2 and 1.25 depending on the bin, much larger than can be reasonably explained by differences in the mass resolution between data and simulation after the calibration using the ψ(2S) data described in Sec. 5. Care is needed in the interpretation of the measured Γ BW and m χ c1 (3872) parameters since m D 0 + m D * 0 − m χ c1 (3872) < Γ BW . The Breit-Wigner parameterisation may not be valid since it neglects the opening of the D 0 D * 0 channel. 7 Flatté model 7

.1 The Flatté lineshape model
The proximity of the χ c1 (3872) mass to the D 0 D * 0 threshold distorts the lineshape from the simple Breit-Wigner form. This has to be taken into account explicitly. The general solution to this problem requires a full understanding of the analytic structure of the coupled-channel scattering amplitude. However, if the relevant threshold is close to the resonance, simplified parametrisations are available and have been used to describe the χ c1 (3872) lineshape [15,16].
In the J/ψ π + π − channel the χ c1 (3872) lineshape as a function of the energy with respect to the D 0 D * 0 threshold, E ≡ m J/ψ π + π − − (m D 0 + m D * 0 ), can be written as where Γ ρ (E) is the contribution of the J/ψ π + π − channel to the width of the χ c1 (3872) state. The complex-valued denominator function, taking into account the D 0 D * 0 and D + D * − two-body thresholds, and the J/ψ π + π − J/ψ π + π − π 0 channels, is given by The Flatté energy parameter, E f , is related to a mass parameter, m 0 , via the relation . The width Γ 0 is introduced in Ref. [16] to represent further open channels, such as radiative decays. The model assumes an isoscalar assignment of the χ c1 (3872) state, using the same effective coupling, g, for both channels. The relative momenta of the decay products in the rest frame of the two-body system, k 1 for D 0 D * 0 and k 2 for the D + D * − channel, are given by where δ = 8.2 MeV is the isospin splitting between the two channels. The reduced masses are given by For m J/ψ π + π − masses below the D 0 D * 0 and D + D * − thresholds these momenta become imaginary and thus their contribution to the denominator will be real. The energy dependence of the J/ψ π + π − and J/ψ π + π − π 0 partial widths is given by [16] Γ The known values for masses m ρ , m ω and widths Γ ρ , Γ ω [31] are used and the lineshapes are approximated with fixed-width Breit-Wigner functions. The partial widths are parameterised by the respective effective couplings f ρ and f ω and the phase space of these decays, where intermediate resonances ρ 0 → π + π − and ω → π + π − π 0 are assumed. The dependence on E is given by the upper boundary of the integrals M (E) = E + (m D 0 + m D * 0 ) − m J/ψ . The momentum of the two-or three-pion system in the rest frame of the χ c1 (3872) is given by The model as specified contains five free parameters: m 0 , g, Γ 0 and the effective couplings f ρ and f ω . In contrast to the Breit-Wigner lineshape, the parameters of the Flatté model cannot be easily interpreted in terms of the mass and width of the state. Instead it is necessary to determine the location of the poles of the amplitude. The analysis proceeds with a fit of the Flatté amplitude to the data, Sec. 7.2, and subsequent search for the poles, Sec. 7.4. The resulting lineshape replaces the Breit-Wigner function described in Sec. 6. It is convolved with the resolution models described in Sec. 4. The Flatté parameters are estimated from a simultaneous unbinned likelihood fit to the J/ψ π + π − mass distribution in the six p π + π − data samples described in Sec. 4. The data points are corrected for the observed shifts of the reconstructed mass of the ψ(2S) in each bin.

Fits of the Flatté lineshape to the data
In order to obtain stable results when using the coupled channel model to describe the J/ψ π + π − mass spectrum, a relation between the effective couplings f ρ and f ω is imposed. This relation requires that the branching fractions of the χ c1 (3872) state to J/ψ ρ 0 and J/ψ ω 0 final states are equal, which is consistent with experimental data [5,14,41], thus eliminating one free parameter in the fit. Furthermore, a Gaussian fit constraint is applied on the ratio of branching fractions The value used here is obtained as the weighted average of the results from the BaBar [5] and Belle [14,41] collaborations, as listed in Ref. [42]. The Flatté model reduces to the Breit-Wigner model as a special case, namely when there is no additional decay channel available near the resonance. However, the R DD * constraint enforces a large coupling to the D 0 D * 0 channel and the lineshape will be different from the Breit-Wigner function in the region of interest. For large couplings to the two-body channel the Flatté parameterisation exhibits a scaling property [43] that prohibits the unique determination of all free parameters on the given data set. Almost identical lineshapes are obtained when the parameters E f , g, f ρ and Γ 0 are scaled appropriately. In particular, it is possible to counterbalance a lower value of E f with a linear increase in the coupling to the DD * channels g. While this is not a true symmetry of the parameterisation -there are subtle differences in the tails of the lineshape -in practice, within the experimental precision this effect leads to strong correlations between the parameters. Figure 3 illustrates the scaling behavior in the data. The black points show the best-fit result for the parameter g evaluated at fixed E f , optimising the remaining parameters at every step. To a good approximation g depends linearly on E f with The red points show the negative log likelihood relative to is minimum value ∆LL for each of these fits, revealing a shallow minimum around m 0 = 3860 MeV. At lower E f values ∆LL raises very slowly, reaching a value of 1 around −270 MeV. Values of E f approaching the D 0 D * 0 threshold are disfavoured, though. In particular good quality fits are obtained only for negative values of E f . A similar phenomenon has been observed in the previous analyses of BaBar and Belle data and is discussed in Ref. [16]. As in those studies, for the remainder of the paper the practical solution of fixing m 0 = 3864.5 MeV, corresponding to E f = −7.2 MeV, is adopted. The remaining model parameters are evaluated with this constraint applied. This procedure has been validated using pseudoexperiments and no significant bias is found. For g and Γ 0 the uncertainties estimated by the fit agree with the spread of the pseudoexperiments. For f ρ an uncertainty which is 10% larger than what is found in the pseudoexperiments is observed and this conservative estimate is reported.
The measured values for g, f ρ and Γ 0 are presented in Table 4. In order to fulfill the constraint on the branching ratios, Eq. (7), the effective coupling, f ω , is found to be 0.01. The systematic uncertainties on the Flatté parameters are summarised in Table 5 and discussed below. The systematic uncertainties introduced by the background and resolution parameterisations are evaluated in the same way as for the Breit-Wigner analysis, using discrete profiling. The impact of the momentum scale uncertainty is investigated by shifting the data points by 66 keV and repeating the fit. Further systematic uncertainties are particular to the Flatté parameterisation. The location of the D 0 D * 0 threshold is known to a precision of 0.11 MeV [31]. Varying the threshold by this amount and repeating the fit leads to an uncertainty on the parameters which is similar to that introduced by the momentum scale. Finally, the D * 0 meson has a finite natural width, for which an upper limit of Γ D * 0 < 2.1 MeV [31] has been measured. However, theoretical predictions estimate Γ D * = 65.5 ± 15.4 keV [37], based on the measured width of the D * + meson. Modified lineshape models taking into account the finite width of the D * 0 are available. In particular, Refs. [37,44] suggest replacing k 1 (E) in Eq. (3) with where E R ≡ m D * 0 − m D 0 − m π 0 . The reduced mass, µ, is calculated as (2m D 0 +m π 0 ) . With this modification there is always a contribution to both the imaginary and real part of the denominator function in Eq. (2). Repeating the fit results in a similar but worse fit [GeV] quality with a log-likelihood difference of 0.1. The width Γ 0 is reduced by 0.2 MeV, which is the smallest systematic uncertainty on this parameter.  Table 6. The mode of the Flatté distribution agrees within uncertainties with the Breit-Wigner solution, discussed in Sec. 6. However, the FWHM of the Flatté model is a factor of five smaller than the Breit-Wigner width. To check the consistency of these seemingly contradictory results, pseudoexperiments generated with the Flatté model and folded with the known resolution function are analysed with the Breit-Wigner model.

Pole search
The amplitude as a function of the energy defined by Eq. (2) can be continued analytically to complex values of the energy E. This continuation is valid up to singularities of the amplitude. There are two types of singularities, which are relevant here: poles and branch points. Poles of the amplitude in the complex energy plane are identified with hadronic states. The pole location is a unique property of the respective state, which is independent of the production process and the observed decay mode. In the absence of nearby thresholds the real part of the pole is located at the mass of the hadron and the imaginary part at half the width of the state. Branch point singularities occur at the threshold of every coupled channel and lead to branch cuts in the Riemann surface on which the amplitude is defined. Each branch cut corresponds to two Riemann sheets. Through Eq. (2) the amplitude will inherit the analytic structure of the square root functions of Eq. (3) that describe the momenta of the decay products in the rest frame of the two-body system. The square root is a two-sheeted function of complex energy. In the following, a convention is used where the two sheets are connected along the negative real axis. An introduction to this subject can be found in Refs. [45][46][47] and a summary is available in Ref. [48]. For the χ c1 (3872) state only the Riemann sheets associated to the D 0 D * 0 channel are important, since all other thresholds are far from the signal region. The following convention is adopted to label the relevant sheets: The fact that the model contains several coupled channels in addition to the D 0 D * 0 channel complicates the analytical structure. The sign in front of the momentum √ −2µ 1 E is the same for sheets I and II and therefore they belong to a single sheet with respect to the D 0 D * 0 channel. The two regions are labelled separately due to the presence of the the J/ψ π + π − , J/ψ π + π − π 0 channels, as well as radiative decays. Those channels have their associated branch points at smaller masses than the signal region. The analysis is performed close to the D 0 D * 0 threshold and points above and below the real axis lie on different sheets with respect to those open channels.
Sheets I and II correspond to a physical sheet with respect to the D 0 D * 0 channel, where the amplitude is evaluated in order to compute the measurable lineshape at real energies E. Sheets III and IV correspond to an unphysical sheet with respect to that channel. Sheet II is analytically connected to sheet IV along the real axis, above the D 0 D * 0 threshold.
In the single-channel case, a bound D 0 D * 0 state would appear below threshold on the real axis and on the physical sheet.
A virtual state would appear as well below threshold on the real axis, but on the unphysical sheet. A resonance would appear on the unphysical sheet in the complex plane [45][46][47]. The presence of inelastic, open channels shifts the pole into the complex plane and turns both a bound state as well as a virtual state into resonances. In the implementation of the amplitude used for the analysis, the branch cut for the D 0 D * 0 channel is taken to go from threshold towards larger energy E, while the branch cuts associated to the open channels Γ(E) are chosen to lie along the negative real axis. The analytic structure around the branch cut associated to the D + D * − threshold is also investigated, but no nearby poles are found on the respective Riemann sheets.
At the best estimate of the Flatté parameters the model exhibits two pole singularities. The first pole appears on sheet II and is located very close to the D 0 D * 0 threshold . The location of this pole with respect to the branch point, obtained using the algorithm described in Ref. [49], is E II = (0.06 − 0.13 i) MeV. Recalling that the imaginary part of the pole position corresponds to half the visible width, it is clear that this pole is responsible for the peaking region of the lineshape. A second pole is found on sheet III. It appears well below the threshold and is also further displaced from the physical axis at E III = (−3.58 − 1.22 i) MeV. Figure 6 shows the analytic structure of the Flatté amplitude in the vicinity of the threshold. The color code corresponds to the phase of the amplitude on sheets I (for Im E > 0) and II (for Im E < 0) in the complex energy plane. The pole on sheet II is visible, as is the discontinuity along the D 0 D * 0 branch cut, which for clarity is also Figure 7: The phase of the Flatté amplitude as obtained from the fit with a finite D * 0 width of Γ D * 0 = 65.5 keV on sheets I (for Im E > −Γ D * 0 /2) and II (for Im E < −Γ D * 0 /2) of the complex energy plane. Since the D * 0 meson is treated as an unstable particle, the D 0 D * 0 branch cut, indicated by the black solid line, is located at Im E = −Γ D * 0 /2. The location of the pole is on the physical sheet with respect to the D 0 D * 0 system. indicated by the black line. The trajectory followed by the pole when taking the limit where the couplings to all channels but D 0 D * 0 are sent to zero is shown in red and discussed below.
As shown in Table 5, taking into account the finite width of the D * 0 has a small effect on the Flatté parameters. However, the analytic structure of the amplitude close to the threshold is changed such that in this case the branch cut is located in the complex plane at Im E = −Γ D * 0 /2. The phase of the amplitude for this case is shown in Fig. 7. The displaced branch cut is highlighted in black. The pole is found at E II = (25 − 140 i) keV in a similar location to the case without taking into account the D * 0 width. In particular, the most likely pole position is on sheet II, the physical sheet with respect to the D 0 D * 0 system. The location of the pole on sheet III is found to be E III = (−3.59 − 1.05 i) MeV, similar to the fit that does not account for D * 0 width.
The uncertainties of the Flatté parameters are propagated to the pole position by generating large sets of pseudoexperiments, sampling from the asymmetric Gaussian uncertainties that describe the statistical and the systematic uncertainties introduced The confidence regions for the location of the poles, corresponding to 68.3%, 95.4% and 99.7% intervals, are shown in Fig. 8. For large values of g the pole on sheet II moves to sheet IV, which is analytically connected to the former along the real axis above threshold. Therefore, sheet II (for Im E < 0) and sheet IV (for Im E > 0) are shown together for this pole. While a pole location on sheet II is preferred by the data, a location on sheet IV is still allowed at the 2σ level. The pole on sheet III is located well below threshold and comparatively deep in the complex plane and is shown on the right-hand side of Fig. 8. For comparison, the location of the confidence region for the first pole on sheets II and IV are also indicated on sheet III.
The positions of both poles depend on the choice of the Flatté mass parameter m 0 . The dependence of the lineshape on m 0 has been explored in the region below threshold and for the results shown in Fig. 3 the corresponding pole positions are evaluated. The location of the pole on sheet II extracted for −17 < E f < 0 MeV are marked by black circles in Fig. 9. For smaller values of m 0 the pole moves closer to the real axis, for values of m 0 approaching the threshold, the pole moves further into the complex plane. For all fits performed the best estimate for the location of the pole is on sheet II. Figure 9 also shows the combined confidence regions, which account for the explored range of E f . For each fit, a sample consisting of 10 5 pseudoexperiments is drawn from the Gaussian distribution described by the covariance matrix of the fit parameters. Only the statistical uncertainties obtained for each fit are used for this study. The resulting samples of pole positions are combined by weighting with their respective likelihood ratios with respect to the best fit. The preferred location of the pole is on sheet II. However, a location of the pole on sheet IV is still allowed at the 2σ level.
The location of the pole on sheet III, in particular its real part, depends strongly on the choice of m 0 . For small values of m 0 the pole moves away from the threshold and has less impact on the lineshape. For m 0 approaching the threshold this pole moves closer to the branch point and closer to the pole on sheet II. Since the asymmetry of the poles with respect to the threshold contains information on the potential molecular nature of the state [50], the values of the pole positions are provided for the most extreme scenario that is still allowed by the data with a likelihood difference of ∆LL = 1, (cf. Fig. 3) at m 0 = 3869.3 MeV. In this case the two poles are found at E II = (0.09 − 0.33 i) MeV and The location of the threshold with respect to the observed location of the peak has a profound impact on the Flatté parameters and therefore on the pole position. The main uncertainties, which affect on which sheet the pole is found, are the knowledge of the momentum scale and the location of the D 0 D * 0 threshold. As shown in Table 5, both effects are of equal importance. Figure 10 shows the statistical uncertainties of the pole on sheet II for the case that the mass scale is shifted up by 66 keV. The pole is moving closer towards the real axis but the preferred location remains on sheet II. A measurement of the lineshape in the D 0 D * 0 channel is needed to further improve the knowledge on the impact of the threshold location.
It is possible to study the behavior of the poles in the limit where only the DD * channels are considered. The trajectory traced by the pole on sheet II when the couplings to the other channels (f ρ , f ω , Γ 0 ) are sent to zero is indicated by the red curve in Fig. 6. The coupling g and the Flatté mass parameter E f are kept fixed while taking this limit. For the best fit solution the pole moves below threshold and reaches the real axis at E = −24 keV staying on the physical sheet with respect to the D 0 D * 0 threshold. This location is consistent with a quasi-bound state in that channel with a binding energy of E b = 24 keV. If the pole lies in the allowed region on sheet IV, taking the same limit also sends the pole onto the real axis below threshold, but on the unphysical sheet with respect to D 0 D * 0 . This situation corresponds to a quasi-virtual state. Both types of solutions are analytically connected along the real axis through the branch cut singularity. Therefore, only upper limits on the binding energy can be set. For the bound state solution and only accounting for statistical uncertainties, the result is E b < 57 keV at 90% confidence level (CL). Including the systematic uncertainties due to the choice of the model this limit becomes E b < 100 keV at 90% CL. Setting the couplings to the other channels to zero causes the pole on sheet III to move to the real axis as well, reaching it at E = −3.51 MeV. The corresponding values extracted at the highest allowed value of m 0 = 3869.3 MeV are E b = 29 keV for the bound state pole and E b = 0.73 MeV for the pole on the unphysical sheet.

Results and discussion
In this paper a large sample of χ c1 (3872) mesons from b-hadron decays collected by LHCb in 2011 and 2012 is exploited to study the lineshape of the χ c1 (3872) meson. Describing the lineshape with a Breit-Wigner function determines the mass splitting between the χ c1 (3832) and ψ(2S) states to be ∆m = 185.588 ± 0.067 ± 0.068 MeV , where the first uncertainty is statistical and the second systematic. Using the known value of the ψ(2S) mass [31] this corresponds to m χ c1 (3872) = 3871.695 ± 0.067 ± 0.068 ± 0.010 MeV , sheet II sheet IV where the third uncertainty is due to the knowledge of the ψ(2S) mass. The result is in good agreement with the current world average [31]. The uncertainty is improved by a factor of two compared to the best previous measurement by the CDF collaboration [4]. The measured value can also be compared to the threshold value, m D 0 + m D * 0 = 3871.70 ± 0.11 MeV. The χ c1 (3872) mass, evaluated from the mean of a fit assuming the Breit-Wigner lineshape, is coincident with the D 0 D * 0 threshold within uncertainties, with δE = 0.01 ± 0.14 MeV. A non-zero Breit-Wigner width of the χ c1 (3872) state is obtained with a value of Γ BW = 1.39 ± 0.24 ± 0.10 MeV.
The values found here for m χ c1 (3872) and Γ BW are in good agreement with a complementary analysis using fully reconstructed B + → χ c1 (3872)K + decays presented in Ref. [51] and combined therein.
Since |δE| < Γ BW , the value of Γ BW needs to be interpreted with caution as coupled channel effects distort the lineshape. To elucidate this, fits using the Flatté parameterization discussed in Refs. [15,16] are performed. The parameters are found to be g = 0.108 ± 0.003 + 0.005 − 0.006 , f ρ = 1.8 ± 0.6 + 0.7 − 0.6 × 10 −3 , Γ 0 = 1.4 ± 0.4 ± 0.6 MeV , with m 0 fixed at 3864.5 MeV. The mode of the Flatté distribution agrees with the mean of the Breit-Wigner lineshape. However, the determined FWHM is much smaller, 0.22 + 0.06 + 0.25 − 0.08 − 0.17 MeV, highlighting the importance of a physically well-motivated lineshape parameterization. The sensitivity of the data to the tails of the mass distribution limits the extent to which the Flatté parameters can be determined, as is expected in the case of a strong coupling of the state to the D 0 D * 0 channel [43]. Values of the parameter E f above −2.0 MeV are excluded at 90% confidence level. The allowed region below threshold is −270 < E f < −2.0 MeV. In this region a linear dependence between the parameters is observed. The slope dg dE f is related to the real part of the scattering length [15] and is measured to be dg dE f = (−15.11 ± 0.16) GeV −1 .
In order to investigate the nature of the χ c1 (3872) state, the analytic structure of the amplitude in the vicinity of the D 0 D * 0 threshold is examined. Using the Flatté amplitude, two poles are found. Both poles appear on unphysical sheets with respect to the J/ψ π + π − channel and formally can be classified as resonances. With respect to the D 0 D * 0 channel, one pole appears on the physical sheet, the other on the unphysical sheet. This configuration, corresponding to a quasi-bound D 0 D * 0 state, is preferred for all scenarios studied in this paper. However, within combined statistical and systematic uncertainties a location of the first pole on the unphysical sheet is still allowed at the 2σ level and a quasi-virtual state assignment for the χ c1 (3872) state cannot be excluded. For the preferred quasi-bound state scenario the 90% CL upper limit of the D 0 D * 0 binding energy E b is found to be 100 keV. The asymmetry of the locations of the two poles, which is found to be substantial, provides information on the composition of the χ c1 (3872) state. In the case of a dominantly molecular nature of a state a single pole close to threshold is expected, while in the case of a compact state there should be two nearby poles [52]. The argument is equivalent to the Weinberg composition criterion [53] in the sense that the asymmetry of the pole location in momentum space determines the relative fractions of molecular and compact components in the χ c1 (3872) wave function [54] |k 2 | − |k 1 | |k 1 | + |k 2 | = 1 − Z .
Here Z is the probability to find a compact component in the wave function. The momentum |k 1 | = 6.8 MeV is obtained by inserting the binding energy of the bound state pole into Eq. (3). The corresponding value for the second pole is |k 2 | = 82 MeV and therefore one obtains Z = 15%. The asymmetry of the poles depends on the choice of m 0 . The asymmetry is reduced as the m 0 parameter approaches the threshold. The largest value for m 0 that is still compatible with the data is 3869.3 MeV. In this case one obtains Z = 33% and therefore the probability to find a compact component in the χ c1 (3872) wave function is less than a third. It should be noted that this argument depends on the extrapolation to the single channel case. For resonances the wave function normalisation used in the Weinberg criterion is not valid and Z has to be replaced by an integral over the spectral density [54]. Nevertheless, the value obtained in this work is in agreement with the results of the analysis of the spectral density using Belle data [41,55], presented in Ref. [16].
The results for the amplitude parameters and in particular the locations of the poles, are systematically limited. In the future, a combined analysis of the χ c1 (3872) → J/ψ π + π − and χ c1 (3872) → D 0 D * 0 channels will make possible improvements to the knowledge on the amplitude parameters.