Study of $\chi_{cJ}$ decaying into $\phi$ $K^*(892)$ $\bar{K}$

Using a data sample of 106 million $\psi(3686)$ events collected with the BESIII detector operated at the BEPCII storage ring, we study for the first time the decay $\chi_{cJ}\to\phi K^{0}_S K^{\pm}\pi^{\mp}$ and $\chi_{cJ}\to\phi K^{+} K^{-}\pi^{0}$ in the E1 radiative transition $\psi(3686)\to\gamma\chi_{cJ}$. The decays are dominated by the three-body decay $\chi_{cJ}\to \phi K^*(892)\bar{K}$. We measure branching fractions for this reaction via the neutral and charged $K^*(892)$ and find them consistent with each other within the expectation of isospin symmetry. In the $K\bar{K}\pi$ invariant mass distribution a structure near the $K^*(892)\bar{K}$ mass threshold is observed, and the corresponding mass and width are measured to be $1412\pm4(\mathrm{stat.})\pm8(\mathrm{sys.}) \mathrm{MeV}/c^2$ and $\Gamma$ = $84\pm12(\mathrm{stat.})\pm40(\mathrm{sys.}) \mathrm{MeV}$, respectively. The observed state favors an assignment to the $h_1(1380)$, considering its possible $J^{PC}$ and comparing its mass, width and decay mode to those reported in the Particle Data Group.

Using a data sample of 106 million ψ(3686) events collected with the BESIII detector operated at the BEPCII storage ring, we study for the first time the decay χcJ → φK 0 S K ± π ∓ and χcJ → φK + K − π 0 in the E1 radiative transition ψ(3686) → γχcJ . The decays are dominated by the threebody decay χcJ → φK * (892)K. We measure branching fractions for this reaction via the neutral and charged K * (892) and find them consistent with each other within the expectation of isospin symmetry. In the KKπ invariant mass distribution a structure near the K * (892)K mass threshold is observed, and the corresponding mass and width are measured to be 1412±4(stat.)±8(sys.) MeV/c 2 and Γ = 84 ± 12(stat.) ± 40(sys.) MeV, respectively. The observed state favors an assignment to the h1(1380), considering its possible J P C and comparing its mass, width and decay mode to those reported in the Particle Data Group.

I. INTRODUCTION
It is well known that the heavy-quark mass provides a natural boundary between the perturbative and nonperturbative regimes. Quarkonium systems are regarded as a unique laboratory to study the interplay between perturbative and nonperturbative effects in Quantum Chromodynamics (QCD). Exclusive quarkonium decays constitute an important basis for investigating these effects. Unlike the J/ψ and ψ(3686), the P-wave charmonia states χ cJ (J = 0, 1, 2) are not directly produced in e + e − collisions, thus are less well understood to date [1]. Obtaining more experimental data on exclusive decays of these χ cJ states is important for a better understanding of their nature and decay mechanisms, as well as for testing QCD based calculations. Exclusive charmonium decays have been investigated widely within QCD. The dominant dynamical mechanism is cc quark annihilation into the minimal number of gluons allowed by symmetries followed by the creation of light quarkantiquark pairs, which form the final state hadrons [2]. The χ cJ states are expected to annihilate into two gluons. Predictions by the color singlet model give smaller decay widths than those determined experimentally [3][4][5], while much better predictions can be obtained if the color octet state is taken into account [6,7]. Since the χ cJ states are produced copiously in the E1 radiative transition of ψ(3686), with branching fractions around 9% [1], the large ψ(3686) data sample taken with the Beijing Spectrometer (BESIII) located at the Beijing Electron-Positron Collider (BEPCII) provides a unique opportunity for detailed studies of χ cJ exclusive decays.
In the quark model, 22 ss sector resonances, collectively called strangeonia, are expected below 2.2 GeV/c 2 . So far only 7 states are widely accepted experimentally, counting the maximally mixed η − η ′ as one ss state [8]. The axial-vector candidate, h 1 (1380), is assigned as the ss partner of the 1 1 P 1 states, h 1 (1170), considering its mass and dominant decay to the K * (892)K final state. Experimentally, h 1 (1380) has been reported in a Partial Wave Analysis (PWA) only by LASS [9] and Crystal Barrel [10]. The nature of this state is still controversial with respect to the predictions made by considering the mixing between SU (3)-singlet and SU (3)-Octet mesons in the 1 3 P 1 and 1 1 P 1 nonets [11] or those made by a relativized quark model [12]. The mass determined by the LASS measurement is significantly smaller than the theory prediction. If the LASS result is confirmed, the h 1 (1380) would seem too light to be the 1 1 P 1 ss member. The Crystal Barrel results are consistent with theory predictions, which means that h 1 (1380) could be a convincing candidate to be the ss partner of the 1 1 P 1 state h 1 (1170). The measurement of the total width of the h 1 (1380) is thought to be complicated by the nearby K * (892)K threshold, where the mass distribution and effective width can not be well described with a traditional Breit-Wigner form. The direct observation of the h 1 (1380) in experiments and the precise measurement of its resonance parameters may shed light on its nature and aid in identifying the ground state axialvector meson nonet in the quark model. Due to conservation of angular momentum and parity, the axialvector strangeonia candidates are not produced in J/ψ (ψ(3686)) radiative decays, but are expected to be produced through the hadronic decay of χ cJ associated with a vector meson φ or in J/ψ (ψ(3686)) decays with a pseudoscalar meson η ′ (η). In this paper, we report the first measurement of the decay χ cJ → φK 0 S K ± π ∓ and χ cJ → φK + K − π 0 in the electric dipole (E1) radiative transition ψ(3686) → γχ cJ . In the following text, if not specified, K * (892)K denotes K * (892) 0K 0 and its isospin-conjugate state K * (892) ± K ∓ , while KKπ denotes both K 0 S K ± π ∓ and K + K − π 0 . The charge conjugated channel is always implied. This analysis is based on a data sample of 1.06 × 10 8 ψ(3686) events collected with the BESIII detector at the BEPCII. Data with an additional integrated luminosity of 44.5 pb −1 [13] at a center-of-mass energy of √ s = 3.65 GeV are used to study continuum contributions.

II. BESIII DETECTOR
The BESIII detector, described in detail in Ref. [14], has an effective geometrical acceptance of 93% of 4π.
It contains a small cell helium-based main drift chamber (MDC) which provides momentum measurements of charged particles; a time-of-flight system (TOF) based on plastic scintillator which helps to identify charged particles; an electromagnetic calorimeter (EMC) made of CsI(Tl) crystals used to measure the energies of photons and provide trigger signals; and a muon system (MUC) made of Resistive Plate Chambers (RPC) located in the iron flux return yoke of the superconducting solenoid. The momentum resolution of the charged particles is 0.5% at 1 GeV/c in a 1 Tesla magnetic field. The energy loss (dE/dx) measurement provided by the MDC has a resolution better than 6% for electrons from Bhabha scattering. The time resolution of the TOF is 80 ps (110 ps) in the barrel (endcaps). The photon energy resolution is 2.5% (5%) at 1 GeV in the barrel (endcaps) of the EMC.
Monte Carlo (MC) simulated events are used to de-termine the detection efficiency, optimize the selection criteria, and study possible backgrounds. A GEANT4based [15] MC simulation software, which includes the geometric and material descriptions of the BESIII detector, the detector response, and digitization models as well as the detector running conditions and performance, is used to generate MC samples. The ψ(3686) resonance is simulated with the kkmc [16] generator, which is an event generator based on precise predictions of the Electroweak Standard Model for the process e + e − → f f + nγ, f = µ, τ, d, u, s, c, b. The beam energy spread and initial state radiation (ISR) are taken into account in the simulation. The subsequent decay processes are generated with EvtGen [18]. Background studies are based on a sample of 10 8 ψ(3686) inclusive decays, generated with the known branching fractions taken from the the Particle Data Group (PDG) [1], or with lundcharm [19] for the unknown decays.

III. EVENT SELECTION
Charged particles are reconstructed from hits in the MDC. Charged tracks are required to be within the acceptance of the MDC, satisfying | cos θ| < 0.93. For each track, the point of closest approach to the interaction point (IP) must be within 1 cm in the plane perpendicular to the beam direction and within ±10 cm along the beam direction. Particle identification (PID) is carried out by combining information from the MDC and TOF. PID probabilities (prob(i)) are calculated for each track according to different particle hypotheses i (i = π, K and p). To be identified as a kaon, a track is required to have prob(K) > prob(π) and prob(K) > prob(p), while pion candidates are required to satisfy prob(π) > prob(K) and prob(π) > prob(p).
Photon candidates are reconstructed from isolated showers in the EMC. Each photon candidate is required to have a minimum energy of 25 MeV in the EMC barrel region (| cos θ| < 0.8) or 50 MeV in the endcap region (0.86 < | cos θ| < 0.92). In order to improve the reconstruction efficiency and the energy resolution, the energy deposited in the nearby TOF counters is included in the photon reconstruction. The timing information from the EMC is used to further suppress electronic noise and energy deposition unrelated to the event of interest.
K 0 S candidates are reconstructed with all combinations of two oppositely charged tracks (without a requirement on the point of closest approach to the IP), assuming both tracks to be pions. A secondary vertex fit is performed for each combination. Each K 0 S candidate is required to have an invariant mass that satisfies |M π + π − − M K 0 S | < 10 MeV/c 2 and a decay length two times larger than its fit error, where M K 0 S is the nominal mass of K 0 S taken from the PDG [1]. If more than one K 0 S is reconstructed within an event, the one with the minimum |M π + π − − M K 0 S | is selected for further analysis.
π 0 candidates are reconstructed from pairs of photons whose invariant mass satisfies [M π 0 − 60] < M γγ < [M π 0 +40] MeV/c 2 , where M π 0 is the nominal mass of π 0 taken from the PDG [1]. An asymmetrical mass window is used for π 0 reconstruction because the photon energy deposited in the EMC has a long tail on the low energy side. A kinematic fit is performed on the selected photon pairs by constraining their invariant mass to the π 0 mass (1C fit). The χ 2 1C from this kinematic fit is required to be less than 25. To remove backgrounds in which the π 0 is falsely reconstructed from a high energy photon paired with a spurious shower, a restriction is imposed on the decay angle | cos θ decay | < 0.95, where θ decay is the polar angle of each photon in the π 0 rest frame with respect to the π 0 direction in the ψ(3686) rest frame. If more than one π 0 is found within an event, the one with the minimum |M γγ − M π 0 | is selected for further analysis.
In the selection of the decay chain ψ(3686) → γχ cJ , χ cJ → φK 0 S K ± π ∓ , φ → K + K − , a candidate event is required to contain a K 0 S candidate, exactly four additional charged tracks with zero net charge, and at least one photon. The four additional charged tracks must be identified as three kaons and one pion according to PID information. In the selection of ψ(3686) → γχ cJ , χ cJ → φK + K − π 0 , φ → K + K − , a candidate event is required to have four charged tracks with zero net charge, one π 0 candidate, and at least one additional photon. The four charged tracks must be identified as two positively charged and two negatively charged kaons, respectively.
To further remove potential backgrounds and to improve the mass resolution, a four-constraint energymomentum conservation kinematic fit (4C fit) is performed.
After the above selection criteria are applied, the decay φ → K + K − is reconstructed using the two oppositely charged kaons whose invariant mass is closest to the nominal φ mass. Fig. 1 shows the K + K − invariant mass of φ candidates for events in which the invariant mass of is between 3.35 and 3.6 GeV/c 2 . A clear φ signal is observed above a very low background. Signal events are extracted by applying a mass window requirement, |M K + K − −M φ | < 15 MeV/c 2 , which corresponds to three times the mass resolution.
The invariant mass distributions of The arrows indicate the φ sideband region.
(K + K − K + K − π 0 ) with the φ mass window requirement are shown in Fig. 2. Significant χ cJ signals are observed with low background. The corresponding scatter plots of the invariant masses of K 0 S π ± versus K ± π ∓ (K + π 0 versus K − π 0 ) are shown in Fig. 3. The dominant processes are the χ cJ → φK * (892)K three body decays.

IV. BACKGROUND ANALYSIS
Since the candidate events are dominated by the three-body decay χ cJ → φK * (892)K, the corresponding branching fractions are measured by imposing a mass window requirement on the Kπ invariant mass, i.e. |M Kπ − M K * (892) | < 100 MeV/c 2 . For convenience, hereafter, an event with χ cJ → φK * (892) ± K ∓ or χ cJ → φK * (892) 0K 0 decay is called a K * (892) ± or K * (892) 0 event. Also, an event satisfying the re- The potential backgrounds for the decay χ cJ → φK * (892)K are studied based on the inclusive MC sample. The following background categories are considered: The first background category contains the nonχ cJ backgrounds, which are dominated by processes such as ψ(3686) → φK * (892)K * (892) and ψ(3686) → γφK * (892)K. MC studies show that these backgrounds do not produce peaks in the χ cJ mass region, and their invariant mass spectrum can be described in the fit with a polynomial function when extracting the χ cJ signal.
The second background category is the non-φ background, which peaks in χ cJ mass region. The main processes of this background are χ cJ → f ′ 0/2 K * (892)K, which have the same final state as that of the signal, where f ′ 0/2 is a scalar or tensor meson, e.g. f 0 (1710), decaying to the K + K − final state. An unbinned maximum likelihood fit is performed to the K + K − invariant mass. The fit result shows that the non-φ background is less than 1.8% in the φ mass window region. When we calculate the number of χ cJ events, a φ sideband as indicated in Fig. 1 is used to estimate the background from non-φ events. The number of non-φ background events is subtracted to obtain the yields of χ cJ signals.
The third background category is composed of non-K * (892) events. From the scatter plots of the invariant mass of K 0 S π ∓ versus K ± π ∓ (K + π 0 versus K − π 0 ) (see Fig. 3), it appears as though the dominant backgrounds FIG. 3. Scatter plot of the invariant mass of (a) K 0 S π ∓ versus K ± π ∓ and (b) K + π 0 versus K − π 0 with the φ mass window requirement. The boxes represent the K * (892) sideband regions described in the text.
are the reactions containing a high mass K * state, e.g. χ cJ → φK * 0 0/2 (1430)K 0 and χ cJ → φK * ± 0/2 (1430)K ∓ . This background is treated similarly to the second background category: these processes have the same final state as the signal and produce peaking backgrounds in the χ cJ mass region (background level 1.6%). Like the signal, χ cJ → φK * (892)K, this background category has horizontal and vertical bands around 1.430 GeV/c 2 , and contaminates in the χ cJ → φK * (892)K, illustrated as boxes 1 and 2 in Fig. 3. The degree of contamination can be evaluated using the K * 2-dimensional (2D) sideband region with |M Kπ − 1425| < 200 MeV/c 2 , illustrated as region 3 in Fig. 3. The K + K − invariant mass for events within the K * 2D sideband region (without the φ mass window requirement) is studied. The K + K − invariant mass peaks around the φ mass region, and the non-φ events in this region are neglected in the fit.
The last background category is the QED background, which is not produced via the ψ(3686) resonance. Backgrounds of this type are evaluated using the data collected at 3.65 GeV and are found to be small and distributed uniformly in the χ cJ mass region. In the fit, the contribution from this type of background is taken into account by the polynomial function for the first back-ground category.
Unbinned maximum likelihood fits are performed to the invariant mass spectra of K + K − K 0 S K ± π ∓ to extract the yields of χ cJ signal in different regions. In the fits, the three χ cJ signals are described by the corresponding MC shapes convoluted with Gaussian functions which represent the difference in resolution between the data and MC. The background is represented by a second order Chebyshev polynomial function. The results of the fits are shown in Fig. 4, and the signal yields are summarized in Table I, where N sig obs (K * ± ) and N sig obs (K * 0 ) are the number of K * (892) ± and K * (892) 0 tagged events, N bkg φ sideband (K * ± ) and N bkg φ sideband (K * 0 ) are those in the φ sideband regions, and N bkg K * sideband in the K * 2-D sideband region.
For the χ cJ → φK + K − π 0 decay mode, only the charged K * (892) ± state is included. In this analysis, the charge conjugate modes are not separated, and the joint branching fraction of χ cJ → φK * (892) ± K ∓ is measured. Figure 5 shows the distributions of the K + K − K + K − π 0 invariant mass for the signal, φ sideband, and K * 2-D sideband regions.
The same fits are performed to the individual distributions. The fit curves are shown in Fig. 5, and the resulting χ cJ signal yields are summarized in Table II, where N sig obs (K * ± ) is the number of events in the signal region, and N bkg φ sideband (K * ± ) and N bkg K * sideband are those in the φ and K * 2-D sideband region. Here, ǫ is the detection efficiency.

VI. BRANCHING FRACTIONS
The branching fractions are calculated using the extracted signal yields that have been corrected according to the detector efficiency and decay branching fractions of intermediate states.
For the χ cJ → φK 0 S K ± π ∓ decay mode, a clear intersection is observed on the scatter plot of the invariant mass of K 0 S π ± versus K ± π ∓ (Fig. 3 (a)). This indicates that some K * (892) ± events are included in the K * (892) 0 tagged events, and vice versa. The relationship between     the observed signal yields and the corresponding branching fractions is given by the following equations: where N ψ(3686) is the number of ψ(3686) events, N sig c/n is the number of K * (892) ±/0 tagged candidates, and Br c/n is the product branching fraction and Br(K 0 S → π + π − ) are taken from the PDG [1], and the branching fractions of χ cJ → φK * (892) ± K ∓ / φK * (892) 0K 0 are to be measured in this analysis. The factor ǫ cc/nn is the detection efficiency for K * (892) ±/0 events to be identified as K * (892) ±/0 tagged events, and ǫ cn/nc is the detection efficiency for K * (892) ±/0 events to be identified as K * (892) ±/0 tagged events. Then the branching fractions are calculated as : Br n = N sig c ǫ cn − N sig n ǫ cc N ψ(3686) (ǫ cn ǫ nc − ǫ nn ǫ cc ) .
The corresponding background subtracted signal yields N sig c/n are calculated as: Here, f φ sideband and f K * sideband are normalization factors; based on the K + K − invariant mass distribution in Fig. 1(a), f φ sideband is taken as 0.37, while f K * sideband is taken to be 1/4. This takes into account the area of the sideband region (box 3 in Fig. 3) relative to that of polluted signal regions (box 1 or 2 in Fig. 3) as well as the existence of both isospin conjugate channels in the sideband regions. The detection efficiencies are evaluated with exclusive signal MC samples, which contain events in the reactions ψ(3686) → γχ cJ , χ cJ → φK * (892) ± K ∓ /φK * (892) 0K 0 with subsequent decay φ → K + K − , K * (892) ± → K 0 π ± or K * (892) 0 → K ± π ∓ ,K 0 /K 0 → K 0 S , K 0 S → π + π − . The ψ(3686) → γχ cJ decay is generated with a 1 + λ cos 2 θ distribution, where θ is the angle between the direction of the radiative photon and that of the positron beam, and λ = 1, -1/3, 1/13 for J = 0, 1, 2, assuming pure E1 transitions. The χ cJ decays are generated with a flat angular distribution. Due to the existence of resonances in the K * (892) ± K ∓ invariant mass, the detection efficiency is evaluated by weighting the signal MC samples by the K 0 S K ± π ∓ invariant mass. This accounts for differences in the distributions between the data and the MC simulation. The detection efficiencies are listed in Table III and the corresponding measured branching fractions Br(χ cJ → φK * (892) ± K ∓ /φK * (892) 0K 0 ) are given in Table VII. For the χ cJ → φK + K − π 0 decay mode, only the charged K * (892) ± intermediate state is included, and the joint branching fractions of the charged conjugate modes are calculated by where Br is the product branching fraction of the other processes in the cascade decay including ψ(3686) → γχ cJ , φ → K + K − , K * (892) ± → K ± π 0 and π 0 → γγ, and all the individual branching fractions are taken from PDG. A factor of 0.25 is determined by the K + K − invariant mass distribution in Fig. 1(b). The K * 2-D background is taken with a weight of 0.5 due to the double area of box 3 in Fig. 3 (b), compared to box 1 or 2. As described above, the detection efficiencies are evaluated using the exclusive MC samples weighted by the K + K − π 0 invariant mass.

VII. SYSTEMATIC UNCERTAINTIES
Several sources of systematic uncertainties are considered in the measurement of the branching fractions. These include differences between the data and the MC simulation for track reconstruction, PID, photon detection, kinematic fitting, π 0 selection, K 0 S reconstruction, the mass window requirement, the fitting process, background estimation, MC modeling, the branching fractions of intermediate state decays, and the luminosity measurement.
(a) Tracking and PID efficiency. The tracking efficiencies for K ± and π ± as functions of transverse momentum have been investigated with the control samples J/ψ → K 0 S K ± π ∓ , K 0 S → π + π − [21] and ψ(3686) → π + π − J/ψ [22]. The uncertainty of the tracking efficiency is 1% for each pion and 1.58% for each kaon. These uncertainties are obtained taking the transverse momentum distributions found in data into account.
The uncertainty due to PID has been studied via the same control samples [22] and is estimated to be 2% for each charged pion and 2.23% for each charged kaon by weighting according to the transverse momentum distributions.
(b) Photon detection efficiency. The uncertainty due to photon detection efficiency is 1% per photon [23]. This is determined from studies of photon detection efficiencies with a large and high purity control sample of J/ψ → ρ 0 π 0 where ρ 0 → π + π − and π 0 → γγ.
(c) 4C-kinematic fit. In the analysis, the track helix parameters (φ 0 , κ, tan λ) are corrected for the MC sample in order to reduce the difference of the 4C kinematic fit χ 2 4C between the data and the MC sample. Here, φ 0 is the azimuthal angle that specifies the pivot with respect to the helix center, κ is the reciprocal of the transverse momentum and tan λ is the slope of the track. The correction factors are obtained from a clean sample of J/ψ → φf 0 (980), φ → K + K − and f 0 (980) → π + π − . An alternative detection efficiency is evaluated with the same MC samples, but without helix parameters corrections. The difference in the efficiencies is taken as the uncertainty of the 4C kinematic fit [25].
(d) π 0 selection. The uncertainty due to the π 0 selection is determined from a high purity control sample of J/ψ → π + π − π 0 . The difference in the π 0 selection efficiency between the data and the MC simulation, 1.0%, is taken as the uncertainty for the π 0 selection [23].
(e) K 0 S reconstruction. The uncertainty for the K 0 S reconstruction efficiency is studied with a control sample of J/ψ → K * (892) ± K ∓ . A conservative value of 3.5% is taken as the systematic uncertainty for K 0 S reconstruction [24]. This uncertainty is dominated by two sources: one is the tracking efficiency of the two pions from the K 0 S decay, and the other is the secondary vertex fit for the two pions and the related selection criteria.
(f ) φ and K * (892) mass window requirement. The uncertainty from the φ mass window requirement is estimated by changing the φ signal and sideband windows by one time the mass resolution of φ, 5MeV, which is determined by a fit to data in Fig. 1 . The difference in the branching fractions is taken as the uncertainty. For the uncertainty related to the K * (892) mass window requirement, the value for the K * (892) width used in the MC simulation, Γ = 50.8 ± 0.18 MeV [1], is changed by twice the uncertainty quoted by PDG, and the difference in the detection efficiency is taken as the systematic uncertainty.
(g) Fitting process. To estimate the uncertainties from the fitting process, three effects are considered. (1) χ cJ signal lineshape. The χ cJ signal lineshapes are described with MC simulated lineshapes convoluted with Gaussian functions. An alternative fit with Breit-Wigner functions convoluted with Gaussian functions for the χ cJ signals is performed, where the Gaussian functions represent the mass resolution. The difference in the production yield returned by the fits is considered as the systematic uncertainty related to the signal lineshape. (2) background lineshape. The non-peaking background is described with a second order Chebyshev polynomial. Alternative fits with different order Chebyshev polynomial functions are performed. The largest difference in the branching fractions is taken as the systematic uncertainty. (3) fitting range. The invariant mass of (h) Peaking backgrounds. The peaking backgrounds without a φ signal are estimated with the events in the φ sideband region [1.10, 1.13] GeV/c 2 . The corresponding uncertainty is studied by changing the φ sideband to the range [1.08, 1.11] GeV/c 2 . The uncertainty of the peaking background with an excited K * is estimated by the difference in the branching fractions with or without the K * background subtraction.
(i) Weighting method in MC. To obtain the detection efficiency, the MC samples are weighted by the K * (892)K invariant mass to compensate for the difference between the data and MC simulation. To get the uncertainty from this weighting method, the weight in each K * (892)K invariant mass bin is randomly changed around its mean value by a Gaussian with a standard deviation given by the statistical uncertainty. The same process is performed ten thousand times, and the standard deviation on the detection efficiencies is taken as the systematic uncertainty related with the weighting method.
(j) Other uncertainties. The uncertainty of the total number of ψ(3686) decays is 0.8%. This uncertainty is determined from a study of inclusive ψ(3686) hadronic decays [20]. The uncertainties due to the branching fractions of intermediate states are taken from the PDG [1]. Table IV and V summarize the systematic uncertainties for the χ cJ → φK 0 S K ± π ∓ and χ cJ → φK + K − π 0 decay modes, respectively. Assuming all of the uncertainties are independent, the total systematic uncertainties are obtained by adding the individual contributions in quadrature. The KKπ invariant mass distributions are studied in order to identify any intermediate states. Figure 6 shows the distributions of the KKπ mass for the candidate events within the χ 0,1,2 signal regions of the data as well as the corresponding phase space MC samples χ cJ → φK * (892)K. The χ c0,1,2 signal regions are defined as [3.365, 3.455] GeV/c 2 , [3.490, 3.530] GeV/c 2 and [3.540, 3.575] GeV/c 2 , respectively. A threshold enhancement, which can not be described with the phase space, is observed in both χ c1 and χ c2 signal regions (Fig. 6 (b), (c), (e), (f)), but is absent in the χ c0 signal region (Fig. 6  (a), (d)).
Assume that the threshold enhancement comes from an intermediate state, called X. Since the X state is produced in the χ c1,2 decay associated with a φ meson, its C parity is must be negative due to the conservation of C parity. The X state is observed in the K * (892)K final state, and in the χ c2 decay process associated with a φ meson, which rules out an assignment of J P = 0 + . If we assume that the X state is a conventional qq state, we can also rule out a J P of 0 − , 2 + etc. Taking into account its mass region, its decay through K * (892)K, and the accessible J P C and comparing with the PDG [1], the most likely assignment for X is the h 1 (1380) (J P C = 1 +− ). In the following, we assume that the K * (892)K threshold enhancement is the h 1 (1380) intermediate state.
Besides the K * (892)K threshold enhancement observed in the χ c1,2 signal regions, a clear structure around 1.7 GeV/c 2 is observed in the χ c0 decay, also evident in the χ c1,2 decay (Fig. 6). Considering the mass region, and its decay through K * (892)K, we conclude that this structure is the φ(1680). The C parity should be minus, and the φ(1680) meets this requirement. There is also a possible φ(1850) contribution visible in Fig. 6.
To extract the resonance parameters of the h 1 (1380), a simultaneous fit is performed to the K 0 S K ± π ∓ and K + K − π 0 invariant mass distributions for the χ c1,2 candidate events (Fig. 6 (b), (c), (e), (f)). Three resonance states, h 1 (1380), φ(1680) and φ(1850) as well as a phase space contribution from χ c1,2 → φK * (892)K are included in the fit without interference. The phase space distribution is described using the shape of the MC samples, smoothed with the KEYS method [27]. The shapes of φ(1680) and φ(1850) are described by relativistic Breit-Wigner functions with constant widths. All the parameters of the resonances are fixed to PDG values [1]. Based on the results of the χ cJ fits, the background from nonχ cJ contributes is about 2.7%. MC studies show that the shape of the KKπ invariant mass of the non-χ cJ contribution is similar to that of the phase space. Therefore, this background is included in the phase space contribution in the fit to the K 0 S K ± π ∓ and K + K − π 0 invariant mass distributions. Since the h 1 (1380) resonance is close to the K * (892)K mass threshold, it is parameterized with a relativistic Breit-Wigner function with a    mass dependent width [28]. A detailed description of the h 1 (1380) lineshape used in the following fits is given in Appendix A. The simultaneous fit to the K 0 S K ± π ∓ and K + K − π 0 invariant mass distributions is performed for the candidate events in the χ c1 and χ c2 mass regions. The fit yields a mass of (1412.4 ± 4.4) MeV/c 2 and a width of (84.5 ± 12.4) MeV for the h 1 (1380) resonance, where the errors are statistical only. In the fit, we take the change of the detection efficiency as a function of the invariant mass into account. The efficiency functions are 0.12 − 0.02 × M K 0 S K ± π ∓ , 0.14 − 0.03 × M K 0 S K ± π ∓ , 0.18 − 0.02 × M K + K − π 0 and 0.18 − 0.03 × M K + K − π 0 for the χ c1 and χ c2 regions, respectively. The product of the relativistic Breit-Wigner functions and the efficiency functions are used to describe the signal. Figure 7 shows the sum of the invariant mass distribution of the K 0 S K ± π ∓ and K + K − π 0 decay modes for the candidate events in the χ c1 and χ c2 mass regions as well as the sum of the corresponding fit curves. The goodness of fit is determined to be χ 2 /n.d.f = 1.09 by projecting all candidate events in 45 bins. The statistical significance of the h 1 (1380) signal is measured to be greater than 10 σ by comparing the likelihood values of the fit with and without the h 1 (1380). The statistical significances of the φ(1680) and φ(1850) signals evaluated by the same method are found to be 4.3σ and 3.2σ, respectively. We consider the interferences between h 1 (1380) and these resonances in the systematic uncertainties.
The different sources of systematic uncertainty for the h 1 (1380) resonance parameters are considered as follows: The KKπ invariant mass distributions, (a) K 0 S K ± π ∓ in the χc0 region; (b) K 0 S K ± π ∓ in the χc1 region; (c) K 0 S K ± π ∓ in the χc2 region; (d) K + K − π 0 in the χc0 region; (e) K + K − π 0 in the χc1 region; (f) K + K − π 0 in the χc2 region. The black dots represent the data, and the histogram shows the phase space MC sample. (b) Background shape. In the fit, the background is described by a smoothed phase space MC sample. An alternative fit is performed using a second-order Chebychev polynomial function to describe the non-resonant. The differences in mass, 4.7 MeV/c 2 , and width, 5.5 MeV, are taken as the systematic uncertainties due to the shape of (f ) φ(1850) resonance. The systematic uncertainty due to the φ(1850) resonance is evaluated in a similar way as for the φ(1680). The largest change in mass and width, 1.1 MeV/c 2 and 3.5 MeV, are taken as the systematic uncertainties related to φ(1850) resonance.
(g) The branching fraction constraint between isospin conjugate decay modes.
In the nominal fit, no constraint is imposed on branching fractions for the the isospin conjugate modes and different final states. The systematic uncertainty related to the branching fraction constraint is estimated by requiring that the branching fraction of h 1 (1380) → K 0 S K ± π ∓ is twice that of h 1 (1380) → K + K − π 0 , which is expected by isospin symmetry. The changes in the mass and width of the h 1 (1380), 0.3 MeV/c 2 and 3.4 MeV, are taken as the systematic uncertainties of this type.
(h) Interference effects. The fit is repeated under the following interference scenarios, (1) the interference between h 1 (1380) and phase space, (2) the interference between h 1 (1380) and φ(1680), (3) the interference between h 1 (1380) and φ(1850), The largest differences to the nominal values of the mass and width, 4.7 MeV/c 2 and 35.3 MeV, are taken as the systematic uncertainties related to interference effects. Table VI shows the systematic uncertainties for the mass and width of h 1 (1380). Assuming that all sources of systematic uncertainty are independent, the total systematic uncertainty is determined to be 8 MeV/c 2 for the mass and 40 MeV for the width.

IX. SUMMARY AND CONCLUSIONS
Using 106 M ψ(3686) events collected with the BE-SIII detector, we present a study of the decays χ cJ → φK 0 S K ± π ∓ and χ cJ → φK + K − π 0 , via the E1 radiative transition ψ(3686) → γχ cJ . The decays are dominated by the three-body reaction χ cJ → φK * (892)K. The branching fractions for this reaction via neutral and charged K * (892) are measured for the first time and are summarized in Tables VII. The measured branching fractions for χ cJ → φK * (892) ± K ∓ in the two different final states are found to be consistent with each other. The ratio of the branching fraction for χ cJ → φK * (892) ± K ∓ to that of χ cJ → φK * (892) 0K 0 is consistent within the expectations of isospin symmetry.
By examining the invariant mass spectrum of KKπ, a significant excess of events above the phase space expectation is observed near the K * (892)K mass threshold in the decays of χ c1,2 , with a significance greater than 10σ. The observed structure has negative C parity, and is expected to be the h 1 (1380) state, considering its mass, width and decay through K * (892)K. A simultaneous fit is performed to the invariant mass distributions of KKπ for the candidate events in the χ c1,2 signal regions. The mass and width of the h 1 (1380) are determined to be 1412 ± 4(stat.) ± 8(sys.) MeV/c 2 and 84 ± 12(stat.) ± 40(sys.) MeV, respectively. This is the first direct observation of the h 1 (1380) in its decay to K * (892)K. Evidence is also found for the decays χ cJ → φφ(1680) and χ cJ → φφ(1850), but with significances less than 5σ. More data and advanced analysis techniques, e.g. PWA, may shed light on the properties of the structures observed in the KKπ invariant mass spectrum. and for S-wave K * (892)K D K * = |BW K * 12 · t 3 + g · BW K * 13 · t 2 |, where BW K * ij are the usual P -wave relativistic Breit-Wigner functions [29] with the mass and width of K * (892) fixed at their PDG values. Here t i are the Zemach vectors describing the spin 1 content, and are the vector part of and p i (i = 1, 2, 3) are the vector parts of the 4-momenta of π, K andK in the KKπ rest system. The interference between the K * (892)K andK * (892)K is fixed by C parity: constructive (g = +1) for C = +1, and destructive (g = −1) for C = −1. In our analysis, C = −1.
To simplify the fitting process, a numerical method is used to extract the integral of the decay matrix element over the available phase space. The integral can be rewritten as : where the term 1 N N i=1 D 2 K * (i) can be extracted using the MC samples of χ cJ → φX (X → KKπ) with a large statistics N, and V is the phase space volume for the decay X → KKπ.
A set of MC samples generated with a different mass m and zero width for the intermediate state X is used to extract 1 N N i=1 D 2 K * (i). The phase space volume V of a three-body decay with system mass m can be derived as [31]: where m 12 is the invariant mass of system of particle 1 and particle 2, and R 2 is the two-body phase space volume, which can directly taken from PDG [1] : The above integral process gives the same result as the integral over the decay area in the Dalitz decay study [32]. The value of 1 N N i=1 D 2 K * (i), the phase space volume V, and the corresponding product, S(m), as function of invariant mass of KKπ system, m, are shown in Fig. 8. Due to the small difference in the mass and width of neu-tral and charged K * (892) and χ c1,2 phase space volume, the calculations are performed for χ c1,2 and the decay mode of K 0 S K ± π ∓ and K + K − π 0 , individually. S(m) represents the mass dependence of the width Γ(m, m 0 ) for the h 1 (1380), which is used to describe the h 1 (1380) line-shape.