Amplitude analysis and branching fraction measurement of $D_{s}^{+} \to K^-K^+\pi^+\pi^+\pi^-$

Using $e^{+}e^{-}$ annihilation data corresponding to a total integrated luminosity of $6.32$ fb$^{-1}$ collected at the center-of-mass energies between 4.178 and 4.226 GeV with the BESIII detector, we perform an amplitude analysis of the decay $D^+_s \to K^-K^+\pi^+\pi^+\pi^-$ and determine the relative fractions and phases of different intermediate processes. Absolute branching fraction of $D^+_s\to K^-K^+\pi^+\pi^+\pi^-$ decay is measured to be ($6.60\pm0.47_{\rm stat.}\pm0.35_{\rm syst.})\times 10^{-3}$. The dominant intermediate process is $D_{s}^{+} \to a_1(1260)^+\phi, \phi\to K^-K^+, a_1(1260)^+\to \rho\pi^+, \rho\to\pi^+\pi^-$, with a branching fraction of $(5.16\pm0.41_{\rm stat.}\pm0.27_{\rm syst.})\times 10^{-3}$.


Introduction
The hadronic decays of D + s mesons are dominated by two-body processes [1], such as D + s → P P , V P , V S, V V , AP and AV , where P , V , S and A denote pseudo-scalar, vector, scalar and axial-vector mesons, respectively. The branching fractions (BFs) of most of these decays can be calculated theoretically [2], even if the non-perturbative contributions, such as final-state interactions, make some of them hard to predict. Therefore, BFs measurements of the D + s two-body decays are important to test the theoretical calculations and can be helpful to understand the decay mechanisms of D + s mesons. Up to now, there are no references to studies focusing on D + s → AV decays. Among them, the process D + s → a 1 (1260) + φ, which is mediated via the diagram in Fig. 1, can be studied in the D + s → K − K + π + π + π − decay. Moreover, experimental study of D + s → K + K − π + π + π − is also helpful to clarify the tension observed in the ratio R(D * ) ≡ B(B → D * τ + ν τ )/B(B → D * + ν ) ( = e, µ), with an average of 0.295 ± 0.011 ± 0.008 provided by the Heavy Flavor Averaging Group, which differs from the Standard Model prediction of 0.258 ± 0.005 by 2.6 standard deviations [3]. This hints a possibile violation of the lepton flavor universality. However, the R(D * ) measurement at LHCb experiment suffers from a large systematic uncertainty due to the limited knowledge of the inclusive D + s → π + π + π − X decay [4,5], where the -1 - decay of B → D * − D + s , D + s → π + π + π − X is main background in the decay chain B 0 → D * − τ + ν τ , τ + → π + π + π − X. A precise measurement of the branching fraction (BF) of D + s → K − K + π + π + π − , which is one of the dominant processes in D + s → π + π + π − X, can provide a useful input to improve the precision of R(D * ).
The E687 [6] and FOCUS [7] experiments reported the fitted yields and BFs of D + s → K − K + π + π + π − relative to D + s → K + K − π + , as shown in Table 1. By performing an analysis of the resonant substructure in the decay D + s → K − K + π + π + π − , the FOCUS experiment observed that the decay proceeds primarily through a quasi-two-body decay involving an a 1 (1260) + resonance. In this paper, we perform the first amplitude analysis, as well as measuring the absolute BF of D + s → K − K + π + π + π − , by using the 6.32 fb −1 data samples collected by the BESIII detector at the center-of-mass energies ( √ s) from 4.178 to 4.226 GeV. More precise measurements and a detailed study of the decay structure are expected. Charge conjugate states are always implied throughout this paper.

Detector and data sets
The BESIII detector [8] records symmetric e + e − collisions provided by the BEPCII storage ring [9], which operates in the center-of-mass energy range from √ s = 2.00 to √ s = 4.95 GeV [10]. The cylindrical core of the BESIII detector covers 93% of the full solid angle and consists of a helium-based multilayer drift chamber (MDC), a plastic scintillator time-of-flight system (TOF), and a CsI(Tl) electromagnetic calorimeter (EMC), which are all enclosed in a superconducting solenoidal magnet providing a 1.0 T magnetic field. The solenoid is supported by an octagonal flux-return yoke with resistive plate counter muon identification modules interleaved with steel. The charged-particle momentum resolution -2 -at 1 GeV/c is 0.5%, and the dE/dx resolution is 6% for electrons from Bhabha scattering. The EMC measures photon energies with a resolution of 2.5% (5%) at 1 GeV in the barrel (end cap) region. The time resolution in the TOF barrel region is 68 ps, while that in the end cap region is 110 ps. The end cap TOF system was upgraded in 2015 using multi-gap resistive plate chamber technology, providing a time resolution of 60 ps [11][12][13].
The data samples used in this analysis contain a total integrated luminosity of 6.32 fb −1 collected at the center-of-mass energies between 4.178 and 4.226 GeV. The integrated luminosity of each data sample is shown in Table 2 [14]. The data sets are organized into three sample groups, 4.178 GeV, 4.189-4.219 GeV, and 4.226 GeV, which were acquired during the same year under consistent running conditions. Inclusive Monte Carlo (MC) samples that are 40 times larger than the data sets are produced at the center-of-mass energies between 4.178 and 4.226 GeV with a geant4based [15] MC package, which includes the geometric description of the BESIII detector and the detector response. These samples are used to determine detection efficiencies and to estimate backgrounds. The samples include the production of open charm processes, the initial-state radiation (ISR) production of vector charmonium(-like) states and the continuum processes incorporated in kkmc [16,17]. The known decay modes are modeled with evtgen [18,19] using the BFs taken from the Particle Data Group (PDG) [1], and the remaining unknown charmonium decays are modeled with lundcharm [20,21]. Final state radiation (FSR) from charged final state particles is incorporated using photos [22]. A phase-space (PHSP) MC sample is produced with the D + s → K − K + π + π + π − generated with a uniform distribution and is used to extract the detection efficiency maps. Initially, the PHSP MC sample is used to calculate the normalization integral used in the determination of the amplitude model parameters in the fit to data. Then, the signal MC sample is regenerated with the D + s meson decaying to K − K + π + π + π − using the amplitude model. It is used to evaluate the fit quality and estimate the systematic uncertainty.

Event selection
The data samples were collected just above the D * ± s D ∓ s threshold. In this energy region, pairs of D * ± s D ∓ s mesons are produced copiously; subsequently, the D * ± s meson predominantly decays to γD ± s with a branching fraction of (93.5 ± 0.7)% [1]. The tag method [23] allows to select clean signal samples, providing the opportunity to perform amplitude analyses and to measure the absolute BFs of the hadronic D + s meson decays. In the tag method, a single-tag (ST) candidate requires only one of the D ± s mesons to be reconstructed via a hadronic decay. A double-tag (DT) candidate has both D + s D − s mesons reconstructed, requiring the D + s meson decaying to the signal mode D + s → K − K + π + π + π − and the D − s meson decaying to the eight tag modes listed in Table 3. The reconstruction of π ± , K ± , K 0 S , π 0 , η, and η particles in the final state is discussed below. Table 3. Requirements on M tag for various tag modes, where the η and η subscripts denote the decay modes used to reconstruct these particles.

Tag mode
Mass window (GeV/c 2 ) All charged tracks reconstructed in the MDC must satisfy the condition |cosθ| < 0.93, where θ is the polar angle with respect to the direction of the positron beam. For charged tracks not originating from K 0 S decays, the distance of closest approach to the interaction point is required to be less than 10 cm along the beam direction and less than 1 cm in the plane perpendicular to the beam. Particle identification (PID) for charged tracks combines the measurements of the dE/dx in the MDC and the time of flight in the TOF to form likelihoods L(h) (h = K, π) for each hadron h hypothesis. Charged kaons and pions are identified by comparing the likelihoods for the two hypotheses, L(K) > L(π) and L(π) > L(K), respectively.
The K 0 S candidates are selected by looping over all pairs of tracks with opposite charges, whose distances to the interaction point along the beam direction are within 20 cm. A primary vertex and a secondary vertex are reconstructed and the decay length between the two vertexes is required to be greater than twice its uncertainty. This requirement is not applied for the D − s → K 0 S K − decay since the combinatorial background is low. Candidate K 0 S particles are required to have the χ 2 of the vertex fit less than 100 and an invariant mass of the π + π − pair (M π + π − ) in the range [0.487, 0.511] GeV/c 2 . To prevent an event being -4 -double counted in the D − s → K 0 S K − and D − s → K − π + π − selections, the value of M π + π − is required to be outside of the mass range [0.487, 0.511] GeV/c 2 for the D − s → K − π + π − decay.
Photon candidates are identified using showers in the EMC. The deposited energy of each shower must be more than 25 MeV in the barrel region (|cos θ| < 0.80) and more than 50 MeV in the end cap region (0.86 < |cos θ| < 0.92). The angle between the position of each shower in the EMC and the closest extrapolated charged track must be greater than 10 degrees to exclude showers that originate from charged tracks. The difference between the EMC time and the event start time is required to be within [0, 700] ns to suppress electronic noise and showers unrelated to the event.
The π 0 (η) candidates are reconstructed through π 0 → γγ (η → γγ) decays, with at least one photon falling in the barrel region. The invariant mass of the photon pair for π 0 and η candidates must be in the ranges [0.115, 0.150] GeV/c 2 and [0.500, 0.570] GeV/c 2 , respectively, which are about three times the mass resolution around their known masses. A kinematic fit that constrains the γγ invariant mass to the π 0 or η known mass [1] is performed to improve the mass resolution and the χ 2 is required to be less than 30. The η candidates are formed from the π + π − η combinations with an invariant mass within a range of [0.946, 0.970] GeV/c 2 .
Eight tag modes are reconstructed and the corresponding mass windows on the tagging D − s mass (M tag ) are listed in Table 3. The signal D + s candidates, whose M rec lies within the mass windows listed in Table 2, are retained for further studies, considering the quantity M rec defined as where E cm is the energy of the initial state calculated from the beam energy, p D − s is the three-momentum of the D − s candidate in the e + e − center-of-mass frame, and m D − s is the D − s known mass [1].

Further selection criteria
The following selection criteria are further applied in order to obtain signal samples with high purity for the amplitude analysis. The selection criteria discussed in this section are not used in the BF measurements. A six-constraint (6C) kinematic fit is performed to the process e + e − → D * ± s D ∓ s → γD + s D − s , assuming D − s decaying to one of the tag modes and D + s decaying to the signal mode (K − K + π + π + π − ) with two hypotheses: the signal D + s comes from a D * + s or the D − s comes from a D * − s . The total four-momentum is constrained to the initial four-momentum of the e + e − system, and the invariant masses of tag D − s and D * ± s candidates are constrained to the corresponding known masses. The best D * ± s D ∓ s combination with the minimum χ 2 6C is selected. Then, a seventh constraint of the signal D + s invariant mass is added to the 6C kinematic fit, in order to ensure that all the events fall into the phase-space (PHSP) boundary. The updated four-momenta obtained from the seven-constraint kinematic fit of the final particles are used to perform the amplitude analysis. The signal region of the D + s invariant mass for the D + s → K − K + π + π + π − decay is defined as [1.955, 1.982] GeV/c 2 . The π + with lower invariant mass of π + π − pair is denoted as π + 1 and the other π + as π + 2 . If the invariant mass of π + 1 π − or π + 2 π − satisfies the selection of K 0 S mesons, these candidates are vetoed. There are background events from D + s → K − K + π + π 0 , π 0 → e + e − γ, in which the e + e − pair is misidentified as a π + π − pair. We use cosθ(π + 1 π − ) > 0.985 to suppress this kind of background, where θ is the opening angle between the momenta of π − and π + 1 . Figure 2 shows the fits to the invariant-mass distributions of the accepted signal D + s candidates (M sig ) for various data samples. The signal is described by a MC-simulated shape convolved with a Gaussian function, and the background is described by a linear function. Finally, a mass window

Fit method
The amplitude analysis of D + s → K − K + π + π + π − decay is performed by using an unbinned maximum likelihood fit. The isobar formulism is used to model the total amplitude. For where p j is the set of the final state particles' four momenta, the index j refers to the different particles in the final states, the indices 1, 2 and 3 correspond to the three intermediate resonances.
In the amplitude, S n (p j ) is the spin factor, The total amplitude M is the coherent sum of the amplitudes of the intermediate processes, M(p j ) = c n A n (p j ), where c n = ρ n e iφn is the corresponding complex coefficient. The magnitude ρ n and phase φ n will be determined in the likelihood fit. The signal probability density function (PDF) f S (p j ) is given by: where (p j ) is the detection efficiency parameterised in terms of the final four-momenta p j and R 5 (p j ) is the PHSP element of five-body decays. In the numerator of Eq. (4.2), (p j ) and R 5 (p j ) terms are independent of the fitted variables, so they are regarded as constant terms in the fit. The normalization integrals are determined by a MC integration: where k MC is the index of the k th MC event of the MC sample, N gen is the number of the generated MC events and N MC is the number of the selected MC events. The M gen (p j ) is the PDF used to generate the MC samples in the MC integration. The computational efficiency of the MC integration is significantly improved by evaluating the normalization integral with signal MC samples, which intrinsically take into account the event selection acceptance and the detection resolution.
The effect from the tracking and PID differences between data and simulation is considered by multiplying the weight of the MC event by a factor γ , which is calculated as: where i refers to tracking or PID, i,data (p j ) and i,MC (p j ) is the tracking or PID efficiency as a function of the momenta of the daughter particles for data and MC, respectively. By weighting each signal MC event with γ , the MC integration is given by: The contribution from the background is subtracted in the likelihood calculation by assigning a negative weight to the background events. The log-likelihood function is written as: 6) where N data is the number of candidate events in data, w bkg k and N bkg are the background weight and the number of simulated background events, respectively.
To combine the data samples taken at various center-of-mass energies, Eq. (4.6) is re-written as: where n denotes the data samples at √ s = 4.178 GeV, 4.189-4.219 GeV, and 4.226 GeV, respectively.

Blatt-Weisskopf barriers
For a decay process a → bc, the Blatt-Weisskopf barriers X L (q) [24] depend on the angular momentum L = 0, 1, 2 and the momentum q of the final-state particle b or c in the rest system of a. They are defined as: with z 0 = q 0 R and z = qR, where R is the effective radius of the intermediate resonances.
The momentum q is given by: where s a , s b and s c refer to the squared invariant masses of particles a, b and c, respectively. The value of q 0 is that of q when s a = m 2 a . The effective radius of barrier R is fixed to be 3.0 GeV −1 for the intermediate resonances and 5.0 GeV −1 for the D + s meson.

Propagator
The intermediate resonances a 1 (1260) and φ are parameterised with the relativistic Breit-Wigner (RBW) formula: where m = E 2 − p 2 , m 0 is the nominal mass of the intermediate resonance, and Γ(m) is given by: -8 -where Γ 0 is the width of the intermediate resonance.
The ρ 0 meson is parameterised with the Gounaris-Sakurai (GS) line shape [25], which is given by: where: 13) and the function h(m) is defined as: with: The normalization condition at P GS (0) fixes the parameter d = f (0) Γ 0 m 0 , which results in: (4.16)

Spin factors
Due to the limited phase space available in the decay, we only consider the states with angular momenta up to 2. As discussed in [26], we define the spin projection operator for a process a → bc, P (S) µ 1 ···µ S ν 1 ···ν S as: The quantities p a , p b , and p c are the momenta of particles a, b, and c, respectively. The covariant tensors are given by: µνµ ν (a)r µ a r ν a . where r a = p b − p c . Eleven kinds of spin factors are listed in Table 4, where the tensor describing the D + s decay is denoted byT and the one of the a decay is denoted byt.

Fit results
The amplitude of the D + s [S] → a 1 (1260) + φ, a 1 (1260) + [S] → ρ 0 π + , φ → K − K + decay is expected to have the largest contribution and it has been chosen as the reference. Thus, its magnitude and phase are fixed to 1.0 and 0.0, respectively, while the other amplitudes are left floating in the amplitude fit. The masses and widths of all the resonances are fixed to the corresponding PDG averages [1]. The background weights are fixed according to the fits shown in Fig. 2.
The calculation of the fit fraction (FF) for each individual amplitude does not involve detector acceptance or resolution effects, and is based on PHSP MC truth information, as it follows: 20) where N gen is the number of PHSP MC generated events,Ã n is either the n th amplitude (Ã n = c n A n ) or the n th component of a coherent sum of amplitudes (Ã n = c n,i A n,i ).
The phases, FFs, and statistical significances for various amplitudes are listed in Table 5. The mass projections of fit results are shown in Fig. 3. The assignments of systematic uncertainties are discussed in the next section.

Systematic uncertainties for amplitude analysis
The systematic uncertainties for the amplitude analysis are summarised in Table 6, with their assignments described below.
i Resonance parameters: the masses and widths of a 1 (1260), φ and ρ are shifted by their corresponding uncertainties [1]; the changes of the phases and FFs are assigned as the associated systematic uncertainties.
ii R values: the systematic uncertainties associated with effective radii of barriers (R values) are estimated by repeating the fit procedure by varying the radii of the intermediate states and D + s meson within 1 GeV −1 .
iii Fit bias: the uncertainty due to the fit procedure is evaluated by studying signal MC samples. An ensemble of 300 signal MC samples are generated according to the nominal results of this analysis. After applying the selection criteria, each of these samples has the same size as the data sample and is used to perform the same amplitude analysis. The pull of each parameter is defined as Out(i)−In(i) , where i denotes the different parameters, In(i) denotes the input value, Out(i) is the value obtained from the fit to a signal MC sample and σ stat(i) is the corresponding statistical uncertainty. For each parameter, 300 pull values are obtained and the deviations of their average from zero are considered as the systematic uncertainty.
iv Background estimation: the background is determined by the inclusive MC samples; the fractions of background events are increased or decreased by one corresponding statistical uncertainties, and the largest differences from the nominal results are considered as the uncertainties.
v Lineshape of the ρ meson: an alternative lineshape parameterization with relativistic Breit-Wigner instead of Gounaris-Sakurai is used, and differences are included inside the uncertainties.
vi Experimental effects: to estimate the systematic uncertainty related to the difference of tracking or PID efficiencies between data and MC simulation, which is γ in Table 6. Systematic uncertainties on the phases and FFs for various amplitudes in units of the corresponding statistical uncertainties. The sources are: (i) the fixed parameters in the amplitudes, (ii) the R values, (iii) fit bias, (iv) background, and (v) shape of the ρ meson.
Amplitude Source Eq. (4.5), the amplitude fit is performed varying the tracking and PID efficiencies according to their uncertainties. However, these uncertainties are estimated to be negligible.

Branching fraction measurement
In addition to the selection criteria for final-state particles described in Sec. 3, for the branching fraction measurement all the pions are subjected to an additional momentum cut p(π) > 100 MeV/c, to remove soft pions from D * + decays. For multiple ST candidates, the candidate with M rec closest to the known mass of D * + s [1] is chosen as the best candidate. Besides the tag modes shown in Table 3 in the amplitude analysis, we add another two tag modes: D − s → π − π − π + and D − s → K − π − π + . The M tag windows are [1.952, 1.982] and [1.953, 1.986] GeV/c 2 , respectively. The yields for various tag modes are listed in Table 7, and they are obtained by fitting the corresponding M tag distributions. To prevent an event being double counted in the D − s → K 0 S K − and D − s → K − π + π − selections, the value of M π + π − is required to be outside of the mass range [0.487, 0.511] GeV/c 2 for the D − s → K − π + π − decay. As an example, the fits to the data sample at √ s = 4.178 GeV are shown in Fig. 4. In the fits, the signal is modeled by a MC-simulated shape convolved with a Gaussian function to take into account the data-MC difference. The background is described by a second-order Chebyshev polynomial. For the tag mode D − s → K 0 S K − , there are some peaking backgrounds coming from D − → K 0 S π − . The shape of this background is taken from the inclusive MC samples and added to the fit leaving its yield floating. For the tag mode D − s → π − η , there is the peaking background coming from D − s → ηπ + π − π − . The shape and yield of this background are taken from the inclusive MC samples and added to the fit.
Once a tag mode is identified, we search for the signal decay D + s → K − K + π + π + π − at the recoiling side. In the case of multiple candidates, the DT candidate with the average mass, (M sig + M tag )/2, closest to the D ± s nominal mass is retained.

Tag mode
To measure the BF, we start from the following equations for one ST mode: 2) where N ST tag is the ST yield for the tag mode, N DT tag,sig is the DT yield, N D + s D − s is the total number of D * ± s D ∓ s pairs produced in the e + e − collisions, B tag and B sig are the BFs of the tag and signal modes, respectively, ST tag is the ST efficiency to reconstruct the tag mode and DT tag,sig is the DT efficiency to reconstruct both the tag and signal decay modes. In the case of more than one tag mode and sample group, where α represents the tag modes in the i th sample group. We isolate B sig by using Eq. (5.1): 4) where N ST α,i and ST α,i are obtained from the data and inclusive MC samples, respectively, and DT α,sig,i is determined with signal MC samples. The decay of D + s → K − K + π + π + π − events is generated according to the results of the amplitude analysis.
The DT yield N DT total is found to be 309 ± 22 from the fit to the M sig distribution of the selected DT candidates, with purity (60.4 ± 2.8)% for the data samples at √ s = 4.178-4.226 GeV. The fit result is shown in Fig. 5, where the signal shape is described by a MC-simulated shape convolved with a Gaussian function and the background shape is described by a linear function. Taking into account the differences in K ± and π ± tracking and PID efficiencies between data and MC simulation, we determine the BF of D + s → K − K + π + π + π − to be (6.60 ± 0.47 stat. ± 0.35 syst. ) × 10 −3 .
-15 - The uncertainty in the total yield of the ST D − s mesons is assigned to be 0.4% by taking into account the background fluctuations in the fit, and by examining the changes of the fit yields when varying the signal and the background shapes. The tracking and PID efficiencies of K ± are studied with e + e − → K + K − K + K − and e + e − → K + K − π + π − (π 0 ) events. The data-MC efficiency ratios of K + (K − ) tracking and PID efficiencies, weighted by the corresponding momentum spectra from signal MC events, are 1.005 ± 0.017 (0.998 ± 0.015) and 0.983±0.003 (0.983±0.003), respectively. After correcting the MC efficiencies for these averaged data-MC differences, the systematic uncertainties of tracking and PID efficiencies per K + (K − ) are assigned as 1.7% (1.5%) and 0.3% (0.3%), respectively. The π ± tracking and PID efficiencies are studied with e + e − → K + K − π + π − events. The data-MC efficiency ratios of the π + (π − ) tracking and PID efficiencies are 0.999 ± 0.005 (0.990 ± 0.005) and 1.004 ± 0.002 (1.004 ± 0.002), respectively, and we assign 0.5% (0.5%) and 0.2% (0.2%) as the systematic uncertainties arising from π ± tracking and PID, respectively.
The systematic uncertainty due to the signal shape is studied by repeating the fit without the convolved Gaussian function. For the background shape of the signal D + s , the MC-simulated shape is used to replace the linear function. The difference of the DT yields is taken as the systematic uncertainty. The uncertainty due to the limited MC statistics is where f i is the tag yield fraction, and i and δ i are the signal efficiency and the corresponding uncertainty of tag mode i, respectively. The uncertainty from the amplitude model is estimated by varying the model parameters based on their error matrix. The distribution of 300 efficiency values resulting from this variation are fitted by a Gaussian function and the fitted resolution divided by the mean value is taken as uncertainty. All of the systematic uncertainties are summarised in Table 8. Adding them in quadrature gives a total systematic uncertainty in the BF measurement of 5.4%.

A All other tested amplitudes
All other tested amplitudes' significances are less than 5 σ, so they are not included in nominal fit and their significances are listed in Table 10.  (980)