Study of the Decays $D_{s}^{+} \rightarrow K_{S}^{0}K^{+}$ and $K_{L}^{0}K^{+}$

Using an $e^{+}e^{-}$ annihilation data sample corresponding to an integrated luminosity of $3.19~\mathrm{fb}^{-1}$ and collected at a center-of-mass energy $\sqrt{s} = 4.178~\mathrm{GeV}$ with the BESIII detector, we measure the absolute branching fractions $\mathcal{B}(D_{s}^{+} \rightarrow K_{S}^{0}K^{+}) = (1.425\pm0.038_{\rm stat.}\pm0.031_{\rm syst.})\%$ and $\mathcal{B}(D_{s}^{+} \rightarrow K_{L}^{0}K^{+}) =(1.485\pm0.039_{\rm stat.}\pm0.046_{\rm syst.})\%$. The branching fraction of $D_{s}^{+} \rightarrow K_{S}^{0}K^{+}$ is compatible with the world average and that of $D_{s}^{+} \rightarrow K_{L}^{0}K^{+}$ is measured for the first time. We present the first measurement of the $K_{S}^{0}$-$K_{L}^{0}$ asymmetry in the decays $D_{s}^{+} \rightarrow K_{S,L}^{0}K^{+}$, and $R(D_{s}^{+} \rightarrow K_{S,L}^{0}K^{+})=\frac{\mathcal{B}(D_{s}^{+} \rightarrow K_{S}^{0}K^{+}) -\mathcal{B}(D_{s}^{+} \rightarrow K_{L}^{0}K^{+})}{\mathcal{B}(D_{s}^{+} \rightarrow K_{S}^{0}K^{+}) +\mathcal{B}(D_{s}^{+} \rightarrow K_{L}^{0}K^{+})}= (-2.1\pm1.9_{\rm stat.}\pm1.6_{\rm syst.})\%$. In addition, we measure the direct $CP$ asymmetries $A_{\rm CP}(D_{s}^{\pm} \rightarrow K_{S}^{0}K^{\pm}) = (0.6\pm2.8_{\rm stat.}\pm0.6_{\rm syst.})\%$ and $A_{\rm CP}(D_{s}^{\pm} \rightarrow K_{L}^{0}K^{\pm}) = (-1.1\pm2.6_{\rm stat.}\pm0.6_{\rm syst.})\%$.


I. INTRODUCTION
Two-body hadronic decays of charmed mesons, D → P 1 P 2 (where P 1,2 denotes a pseudoscalar meson), serve as an ideal environment to improve our understanding of the weak and strong interactions because of their relatively simple topology [1,2]. Charmed-meson decays into hadronic final states that contain a neutral kaon are particularly attractive. Bigi and Yamamoto [3] first pointed out that the interference of the decay amplitudes of the Cabibbo-favored (CF) transition D →K 0 π and the doubly-Cabibbo-suppressed (DCS) transition D → K 0 π can result in a measurable K 0 S -K 0 L asymmetry A similar asymmetry can be defined in D + s decays by replacing π with K. Additionally, as pointed out in Ref. [4], the interference between CF and DCS amplitudes can also lead to a new CP violation effect, which is estimated to be of an order of 10 −3 . The measurement of K 0 S -K 0 L asymmetries and CP asymmetries in charmed-meson decays with neutral kaon provides insight into the DCS process, as well as information to explore D 0 -D 0 mixing, CP violation and SU(3) flavor-symmetry breaking effects in the charm sector [5,6].
On the theory side, different phenomenological models give predictions for the K 0 S -K 0 L asymmetries: the topological-diagrammatic approach [2] under the SU(3) flavor symmetry (DIAG) or incorporating the SU(3) breaking effects (SU(3) FB ) [7,8,10], the QCD factorization approach (QCDF) [9], and the factorization-assisted topological-amplitude (FAT) [11]. The predicted K 0 S -K 0 L asymmetries in charmed-meson decays from these different approaches, as well as the measured values reported by the CLEO Collaboration [12] are summarized in Table I. Considering the large range of values predicted for the K 0 S -K 0 L asymmetries, their measurements provide a crucial constraint upon models of the dynamics of charmed meson decays.
Experimentally, D +(0) decays have been studied intensively in the past two decades [13]. However, existing measurements of charmed-strange meson decays suffer from poor precision due to the limited size of available data samples and a relatively small production cross section in e + e − annihilation [14]. The most recent measurement of B(D + s → K 0 S K + ) = (1.52 ± 0.05 stat. ± 0.03 syst. )% was reported by the CLEO Collaboration [15]; the result was obtained using a global fit to multiple decay modes reconstructed in an e + e − annihilation sample corresponding to an integrated luminosity of 586 pb −1 at a center-of-mass energy √ s = 4.17 GeV. The Belle Collaboration reported a measurement of the branching fraction B(D + s →K 0 K + ) (ignoring the contribution from K 0 K) [16] using a data sample corresponding to an integrated luminosity of 913 fb −1 collected at √ s around the Υ(4S) and Υ(5S) resonances. Neither B(D + s → K 0 L K + ) nor the corresponding K 0 S -K 0 L asymmetry have been measured yet.
In this paper, measurements of the absolute branching fractions for the decays D + s → K 0 S K + and D + s → K 0 L K + , the K 0 S -K 0 L asymmetry, and the corresponding CP asymmetries are performed using a sample of e + e − annihilation data collected at √ s = 4.178 GeV with 1.2 ± 0. 6 -the BESIII detector at the BEPCII. The data sample corresponds to an integrated luminosity of 3.19 fb −1 . Throughout the paper, charge conjugation modes are implicitly implied, unless otherwise noted.

II. BESIII DETECTOR AND MONTE CARLO SIMULATION
The BESIII detector is a magnetic spectrometer that operates at the BEPCII e + e − collider [17]. The detector has a cylindrical geometry that covers 93% of the 4π solid angle and consists of several subdetectors. A main drift chamber (MDC) with 43 layers surrounding the beam pipe measures momenta and specific ionization of charged particles. Plastic scintillator time of flight counters (TOF), located outside of the MDC, provide charged-particle identification information, and an electromagnetic calorimeter (EMC), consisting of 6240 CsI(Tl) crystals, detects electromagnetic showers. These subdetectors are immersed in a magnetic field of 1 T, produced by a superconducting solenoid, and are surrounded by a multi-layered resistive-plate chamber (RPC) system (MUC) interleaved in the steel flux return of the solenoid, providing muon identification. In 2015, BESIII was upgraded by replacing the two endcap TOF systems with multi-gap RPCs, which achieve a time resolution of 60 ps [18]. A detailed description of the BESIII detector is presented in Ref. [19].
The performance of the BESIII detector is evaluated using a geant4-based [20] Monte Carlo (MC) program that includes a description of the detector geometry as well as simulating its response. In the MC simulation, the production of open charm processes directly produced via e + e − annihilation are modeled with the generator conexc [21], which includes the effects of the beam energy spread and initial-state radiation (ISR). The ISR production of vector charmonium states (ψ(3770), ψ(3686) and J/ψ) and the continuum processes (qq, q = u, d, s) are incorporated in kkmc [22]. The known decay modes are generated using evtgen [23], which assumes the branching fractions reported in Ref. [13]; the fraction of unmeasured decays of charmonium states is generated with lundcharm [24]. The final-state radiation (FSR) from charged tracks is simulated by the photos package [25]. A generic MC sample with equivalent luminosity 35 times that of data is generated to study the background. It contains open charm processes, the ISR return to charmonium states at lower mass, and continuum processes (quantum electrodynamics and qq). The signal MC samples of 5.2 million e + e − → D * ± s D ∓ s events are produced; in these samples the D * ± s decays into γ/π 0 /e + e − D ± s , while one D s decays into a specific mode in Table II and the other into the final states of interest K 0 S K ± or K 0 L K ± . The signal MC samples are used to determine the distributions of kinematic variables and estimate the detection efficiencies.

III. DATA ANALYSIS
The cross section to produce e + e − → D * ± s D ∓ s events at √ s = 4.178 GeV is (889 ± 59 stat. ± 47 syst. ) pb, which is one order of magnitude larger than that to produce e + e − → D + s D − s events [14]. Furthermore, the decay branching fraction B(D * + s → γD + s ) is (93.5 ± 0.7)% [13]. Therefore, in the data sample used, D + s candidates arise mainly from the process e + e − → D * ± s D ∓ s → γD + s D − s , along with small fractions from the processes e + e − → D * ± s D ∓ s → π 0 D + s D − s and e + e − → D + s D − s . The outline of the reconstruction is described first, with all details given later in this section.
In this analysis, a sample of D − s mesons is reconstructed first, which are referred to as "single tag (ST)" candidates. The ST candidates are reconstructed in 13 hadronic decay modes that are listed in Table II Here, π 0 and η candidates are reconstructed from a pair of photon candidates, K 0 S candidates are formed from π + π − pairs, and ρ ±(0) candidates are reconstructed from π ± π 0(∓) pairs, unless otherwise indicated by a subscript.
In the sample of events with an ST candidates, the process D + s → K 0 S K + is reconstructed by selecting a charged kaon and a K 0 S candidates from those not used to reconstruct the ST candidates, which is referred as "double tag (DT)". To reconstruct the D + s → K 0 L K + decay, the photon from the decay D * ± s → γD ± s and the charged kaon from D + s decay are selected to determine the missing-mass-squared where P e + e − is the four-momentum of the e + e − initial 141285±631 42.15±0.03 13.58±0.07 16.33±0.10 K − π + π − 18051±575 48.84±0.26 16.35±0.08 19.73±0.12 π + π − π − 40573±964 56.05±0.18 18.47±0.08 22.55±0.12 47±0.06 9.70±0.08 state and P i (i = D − s , γ, K + ) is the four-momentum of the corresponding particle.
Ignoring the small contribution from the process is the total number of e + e − → D * ± s D ∓ s events in the data sample, B i tag is the branching fraction for the i th ST decay mode, B sig is the branching fraction of the signal decay; ǫ i ST and ǫ i DT are the ST and DT detection efficiencies, respectively, which are evaluated from the signal MC samples corresponding to the i th tag mode. The value of ǫ i DT includes the branching fraction B(K 0 S → π + π − ) of signal side in the analysis of D + s → K 0 S K + . The factors of two in Eqs. (3) and (4) are the result of including chargeconjugated modes in the analysis. We combine Eqs. (3) and (4) for each of the 13 tag modes to obtain where N tot DT = i N i DT is the total number of DT events.

A. Selection of ST events
Good charged tracks, except for the daughter tracks of K 0 S candidates, are selected by requiring the track trajectory approaches the interaction point (IP) within ±10 cm along the beam direction and within 1 cm in the plane perpendicular to the beam direction. In addition, the polar angle θ between the direction of the charged track and the beam direction must be within the detector acceptance by requiring | cos θ| < 0.93. Charged particle identification (PID) is performed by combining the ionization-energy loss (dE/dx) measured by the MDC and the time-of-flight measured by the TOF system. Each charged track is characterized by the PID likelihood for the pion and kaon hypotheses, which are L(π) and Good photon candidates are selected from isolated electromagnetic showers which have a minimum energy of 25 MeV in the EMC barrel region (| cos θ| < 0.8) or 50 MeV in the EMC endcap region (0.86 < | cos θ| < 0.92). To reduce the number of photon candidates that result from noise and beam backgrounds, the time of the shower measured by the EMC is required to be less than 700 ns after the beam collision. The opening angle between a photon and the closest charged track is required to be greater than 10 • , which is used to remove electrons, hadronic showers and photons from FSR. π 0 and η → γγ candidates are reconstructed from pairs of photon candidates that have an invariant mass within the intervals (0.115, 0.150) and (0.50, 0.57) GeV/c 2 , respectively. To improve the momentum resolution, a kinematic fit is performed, constraining the γγ invariant mass to its nominal value [13]; the χ 2 of the fit is required to be less than 20 to reject combinatorial background. η → π + π − π 0 candidates are selected by requiring the corresponding invariant mass to be within the interval (0.534, 0.560) GeV/c 2 .
In order to improve the efficiency of the K 0 S selection, K 0 S candidates are reconstructed from tracks assumed to be pions without PID, and the daughter tracks are required to have a trajectory that approaches the IP to within ±20 cm along the beam direction and | cos θ| < 0.93. The K 0 S candidates are formed by performing a vertex-constrained fit to all oppositely-charged track pairs. To suppress combinatorial background, the χ 2 of the vertex fit is required to be less than 200 and a secondary vertex fit is performed to ensure that the K 0 S candidate originates from the IP. The flight length L, defined as the distance between the common vertex of the π + π − pair and the IP in the plane perpendicular to beam direction, is obtained in the secondary vertex fit, and is required to satisfy L > 2σ L , where σ L is the estimated uncertainty on L; this criterion removes combinatorial background formed from tracks originating from the IP. The four-momenta after the secondary vertex fit are used in the subsequent analysis. The K 0 S candidate is required to have a mass within the interval (0.487, 0.511) GeV/c 2 .
To suppress the background with D * decay D * → πD, the momentum of charged and neutral pions is required to be greater than 100 MeV/c. For K − π + π − ST candidates, the invariant mass of the π + π − pair is required to be outside the interval (0.478, 0.518) GeV/c 2 to remove D − s → K 0 S K − decays. The ST D − s candidates are reconstructed via all the possible selected particles combinations.
The invariant mass of the system recoiling against the selected D − s is defined as where p is the momentum of the ST D − s candidate in e + e − CM frame, and M Ds is the nominal mass of the D s meson [13]. M rec is required to be within the interval (2.05, 2.18) GeV/c 2 . For a specific ST mode, if there are multiple combinations satisfying the selection criteria, only the candidate with the minimum value of |M rec − M D * s | is retained for further analysis. These requirements also accept the events in which the ST D s comes from the decay of the primary D * s . To determine the ST yield, a binned maximum likelihood fit to the distribution of the D − s invariant mass M tag is performed for each tag mode; the distributions and fit results are shown in Fig. 1. In the fit, the probability density function (PDF) that describes the signal is the shape of the signal MC distribution, taken as a smoothed histogram and convolved with a Gaussian function to account for any resolution difference between data and MC simulation. The background is described by a second or third-order Chebychev polynomial function. The ST yields determined by the fits, along with the corresponding ǫ i ST estimated from the generic MC sample, are summarized in Table II.
The signal decay D + s → K 0 S K + is reconstructed recoiling against the selected ST D − s candidate. We select a D + s → K 0 S K + candidate if there is only one K 0 S candidate and one good track, which is identified as a kaon and has positive charge, recoiling against the ST D − s candidate; K + and K 0 S candidates are selected by applying the selection criteria described in Sec. III A. In addition, to suppress combinatorial backgrounds, we reject events in which there are additional charged tracks that satisfy | cos θ| < 0.93 and approach the IP along the beam direction within ±20 cm.
To determine the DT signal yield, a two-dimensional (2D) unbinned maximum likelihood fit is performed on the invariant mass of the K 0 S and K + (M K 0 S K + ) versus M tag distribution of selected events, which is summed over the 13 ST modes, as shown in Fig. 2. In the fit, the total PDF is described by summing over the individual PDFs for the following signal and background components, where x represents M K 0 S K + , and y stands for M tag .
• Signal: F sig (x, y) ⊗ G(x; µ x , σ x ) ⊗ G(y; µ y , σ y ) F sig (x, y) is a 2D function derived from the signal MC distribution by using a smoothed 2D histogram; G(x; µ x , σ x ) and G(y; µ y , σ y ) are Gaussian functions that compensate for any resolution difference between data and MC simulation for the variables M K 0 S K + and M tag , respectively. In the 2D fit, the parameters of G(x; µ x , σ x ) and G(y; µ y , σ y ) are fixed to the values determined by fitting the corresponding one-dimensional (1D) distributions.
• BKGI: F BKGI (x, y) ⊗ G(y; µ y , σ y ) This PDF describes the background composed of a correctly reconstructed ST D − s recoiling against combinatorial background, which are distributed in the horizontal band in Fig. 2. F BKGI (x, y) is derived from the distribution of this type of background in the generic MC sample by using a kernel density estimation method (KEYS) [26]. The resolution function G(y; µ y , σ y ) is the same as that in the Signal PDF.
• BKGII: F BKGII (x, y) ⊗ G(x; µ x , σ x ) This PDF describes the background composed of an incorrectly reconstructed ST D − s recoiling against a correctly reconstructed signal candidate, which are distributed in the vertical band in Fig. 2. F BKGII (x, y) is derived from the distribution of this type of background in the generic MC sample by using KEYS. The resolution function G(x; µ x , σ x ) is the same as that in the Signal PDF.
• BKGIII: P BKGIII (x) × P BKGIII (y) This PDF describes the combinatorial background composed of events in which neither the ST D − s nor signal D + s candidate is correctly reconstructed. These background events do not have any peaking components in either variable. Therefore, BKGIII events are described by two independent secondorder polynomials, P BKGIII (x) and P BKGIII (y), with their parameters determined by the fit to data.
The 2D fit gives a signal yield of 1782 ± 47, where the uncertainty is statistical. The M K 0 S K + and M tag distributions for the data, with the projections of the fit results superimposed, are shown in Fig. 3.   (5), the branching fraction is determined to be B(D + s → K 0 S K + ) = (1.425 ± 0.038 stat. )%.

C. Branching fraction measurement of
The D + s → K 0 L K + candidates are reconstructed by requiring the event has only one good track recoiling against the ST D − s candidate; the charged track is required to be identified as a kaon and have opposite charge compared with ST D − s . The K + is selected with the criteria described in Sec. III A. We further suppress combinatorial backgrounds by requiring no additional charged tracks that statisfy the requirements described in Sec. III B.
In this analysis, the ST and signal candidates are assumed to originate from the decay chain e + e − → D * ± s D ∓ s → γD + s D − s , with one D − s decaying into any of ST modes, and the other decaying into K 0 L K + . We re-construct the K 0 L candidate using a kinematic fit that applies constraints arising from the masses of the ST D − s candidate, the signal D + s candidate, the intermediate state D * ± s , and the initial four-momenta of the event. In the kinematic fit, the K 0 L signal candidate is treated as a missing particle whose four-momentum is determined by the fit. The fit is performed to select γ candidate from the decay D * ± s → γD ± s under two different hypotheses that constrain either the invariant mass of the selected γ and signal D + s or the selected γ and the ST D − s to the nominal mass of the D * − s meson; the hypothesis that results in the minimum value of χ 2 is assumed to be the correct topology. If there are multiple photon candidates, which are not used to reconstruct the ST candidate, the fit is repeated for each candidate and the photon that results in the minimum value of the χ 2 is retained for further analysis. For each event, the four-momentum of the missing particle assumed in the kinematic fit is used to determine the M M 2 of the K 0 L candidate. In order to reduce combinatorial background, χ 2 < 40 is required. To further suppress background with multiple photons, we reject those events with additional photons which have an energy larger than 250 MeV and an opening angle with respect to the direction of missing particle greater than 15 • .
To determine the signal yield, an unbinned maximum likelihood fit is performed on the M M 2 distribution of selected events from all 13 ST modes combined, as shown in Fig. 4. In the fit, three components are included: signal, peaking, and non-peaking backgrounds. The PDFs of these components are described below, where x represents M M 2 .
is derived from the signal MC distribution as a smoothed histogram, and G(x; µ ′ x , σ ′ x ) is a Gaussian function that accounts for any resolution difference between data and MC simulation. The value of σ ′ x is fixed in the data fit to the value obtained from a fit to the M M 2 distribution obtained from a D + s → K 0 S K + control sample where the K 0 S is ignored in the reconstruction.
• Peaking background: MC simulated events by using a smoothed histogram. These events form a peaking background if the K 0 S or η is not reconstructed. Here, G(x; µ ′ x , σ ′ x ) is the Gaussian resolution function, whose parameters are the same as those used in the signal PDF. The expected yields of D + s → K 0 S K + and D + s → ηK + are fixed to 263 and 57, respectively. The expected peaking background yields are estimated by using the equation ǫ MC DT are the detection efficiencies of the nominal analysis and the DT method for each mode, respectively; these are estimated from MC simulation samples. The uncertainties of estimated event numbers for D + s → K 0 S K + and D + s → ηK + are 19 and 12, which will be used in systematic uncertainty study.
• Non-peaking background: P (x) P (x) is a function to describe the combinatorial background, which is not expected to peak in the M M 2 distribution. P (x) is a second-order poly- nomial function whose parameters are determined from the fit to data.
The fit to the M M 2 distribution is shown in Fig. 4. The signal yield determined by the fit is 2349 ± 61 events, where the uncertainty is statistical. Using Eq. (5), the branching fraction is calculated to be B(D + s → K 0 L K + ) = (1.485 ± 0.039 stat. )%, where the DT detection efficiencies ǫ K 0 L MM 2 used are summarized in Table II; the values of ǫ K 0 L MM 2 are estimated from signal MC samples.

D. Asymmetry measurement
By using the measured branching fractions and Eq. (1) the K 0 S -K 0 L asymmetry is determined to be To determine the direct CP violation, we also measure the branching fractions for the D + s and D − s decays separately, using the same methodology as the combined branching fraction measurement. The direct CP asymmetriy is defined as which leads to the measurements (10) for the two signal modes.

IV. SYSTEMATIC UNCERTAINTY
For the absolute branching fractions, which are determined according to Eq. (5), the systematic uncertainties are associated with N i ST , N tot DT , and the corresponding ratio of detection efficiencies (ǫ i DT /ǫ i ST ). One of the advantages of the DT method is that most of the systematic uncertainties associated with selection criteria for the ST side reconstruction cancel. However, there is some residual uncertainty due to the different decay topologies between DT and ST events; this is referred to as "tag-side bias", and its effect is considered as one of the systematic uncertainties. For the R(D + s ) and A CP measurements, the systematic uncertainties are calculated by propagating corresponding branching fraction uncertainties from different sources taking into account that some of the uncertainties cancel due to the fact that these observables are ratios as defined in Eqs. (1) and (8). Table III summarizes the relative uncertainties on the absolute branching fraction and the absolute uncertainties for the asymmetries. The total systematic uncertainties are caculated as the sum in quadrature of individual contributions by assuming the sources are independent of one another.
The K + and K − tracking efficiencies are studied using a control sample of e + e − → K + K − π + π − events; the efficiency is calculated as a function of the transverse momentum of the particles. The average efficiency difference between data and MC is computed to be 0.5% by weighting the efficiency difference found in the control sample according to the transverse momentum of kaon in signal MC samples. This is assigned as the systematic uncertainty from this source.
The K + and K − PID efficiencies are studied using a control sample of D + s → K + K − π + , D 0 → K − π + and D 0 → K − π − π + π + events; the efficiency is calculated as a function of the momentum of the particle. The average efficiency difference between data and MC is computed to be 0.5% by weighting the efficiency difference found in the control sample according to the momentum of kaon in signal MC samples, and this is assigned as the systematic uncertainty from this source.
The K 0 S reconstruction efficiency has been studied using control samples of J/ψ → K * (892) ∓ K ± and J/ψ → φK 0 S K ± π ∓ in different momentum intervals [27]. The efficiency difference between data and MC is computed to be 1.5%, which is assigned as the systematic uncertainty from this source.  The systematic uncertainty associated with the photon selection efficiency and the kinematic fit in the study of D + s → K 0 L K + is estimated from the control sample D + s → K + K − π + . The same kinematic fit as that used on the data is performed by assuming the K − π + system is missing. The efficiency difference found between data and MC simulation, 2.0%, is taken as the systematic uncertainty.
The systematic uncertainties associated with the requirements on the energy of additional photons and the number of extra charged tracks are estimated from the control sample D + s → K + K − π + . The efficiency differences between data and MC simulation for these two requirements are both 0.6%, which are assigned as the systematic uncertainties from these sources.
The uncertainty related to the limited sizes of MC samples is 0.3% for both D + s → K 0 S K + and D + s → K 0 L K + . The uncertainties associated with ST, DT, and M M 2 fits are studied by changing the signal and background PDFs, as well as the fit interval; each change is applied separately. Furthermore, in the M M 2 fit, the effect of the assumed peaking background yields is estimated by changing the fixed numbers of events by ±1σ. The systematic uncertainties related to the ST, DT, and M M 2 fit procedure are 0.9%, 0.8% and 1.5%, respectively; these are the sums in quadrature of the relative changes of signal yield that result from each individual change to the fit procedure.
As discussed previously, the selected ST D − s sample is dominated by the process e + e − → D * ± s D ∓ s → γD + s D − s , but there is small contribution from the processes e + e − → D * ± s D ∓ s → π 0 D + s D − s and e + e − → D + s D − s . In the analysis of D + s → K 0 S K + , detailed MC studies indicate that ǫ i DT /ǫ i ST is almost the same for the three processes, since distributions of the kinematic variables are similar and no kinematic fit is performed in the DT selection. Thus, the effect from includ-ing e + e − → D * ± s D ∓ s → π 0 D + s D − s and e + e − → D + s D − s processes is negligible in the absolute branching fraction measurement. In the analysis of D + s → K 0 L K + , the kinematic fit is performed under the hypothesis that the event is e + e − → D * ± s D ∓ s → γD + s D − s , and the MC studies indicate that the contribution of e + e − → D * ± s D ∓ s → π 0 D + s D − s and e + e − → D + s D − s in signal events can be neglected. Thus, the uncertainty of branching fraction B(D * + s → γD + s ) [13] used in the signal MC simulation must be taken as a source of systematic uncertainty. The systematic uncertainty from excluding the process e + e − → D + s D − s is 0.4%, which is the fraction of the ST yields that comes from the process e + e − → D + s D − s ; this fraction is estimated from the MC simulation.
The tag-side bias uncertainty is defined as the uncanceled uncertainty in tag side due to different track multiplicities in generic and signal MC samples. By studying the differences of tracking and PID efficiencies between data and MC in different multiplicities, the tagside bias systematic uncertainties are estimated to be 0.3% for D + s → K 0 S K + and 0.5% for D + s → K 0 L K +

V. SUMMARY AND DISCUSSION
In summary, by using an e + e − collision data sample at √ s = 4.178 GeV, corresponding to an integrated luminosity of 3.19 fb −1 , the absolute branching fractions are measured to be B(D + s → K 0 S K + ) = (1.425 ± 0.038 stat. ± 0.031 syst. )% and B(D + s → K 0 L K + ) = (1.485 ± 0.039 stat. ± 0.046 syst. )%, the former is one standard deviation lower than the world average value [13], and the latter is measured for the first time. The K 0 S -K 0 L asymmetry in D + s decay is measured for the first time as R(D + s → K 0 S,L K + ) = (−2.1 ± 1.9 stat. ± 1.6 syst. )%. This mearsurement is compatible with theoretical predictions listed in Table I. Direct CP asymmetries of the two processes are obtained to be A CP (D ± s → K 0 S K ± ) = (0.6 ± 2.8 stat. ± 0.6 syst. )% and A CP (D ± s → K 0 L K ± ) = (−1.1 ± 2.6 stat. ± 0.6 syst. )%. No significant asymmetries are observed and the uncertainties are statistically dominant.

VI. ACKNOWLEDGMENTS
The BESIII collaboration thanks the staff of BEPCII and the IHEP computing center and the supercomputing center of USTC for their strong support. This work