\section{Introduction}

The $D_{s}^{+}$ meson is the ground state of charmed-strange mesons, and
precise measurements of its semileptonic decays allow for crucial
tests of Standard Model predictions of flavor-changing
interactions. This article presents a new measurement of the
$D_{s}^{+}$-meson inclusive semielectronic branching fraction and positron
momentum spectrum. (Here and throughout this article charge conjugate
modes are implied.)  A previous measurement by the CLEO-c
experiment reported the first results for these quantities,
including the measurement of
$\BF\left(\semilep\right)=\left(6.52\pm0.39\;(\text{stat.})\pm
0.15\;(\text{syst.})\right)\%$~\cite{ref:CLEOMeas}.  Measuring this
branching fraction with improved precision contributes to a
comprehensive understanding of $D_s^+$ decays and is an important
component of the overall experimental and theoretical heavy-flavor
physics program.

\renewcommand{\arraystretch}{1.25}
\begin{table*}[hbt]
\centering{
\begin{tabular}{c|cc}
Mode                                         &  \multicolumn{2}{c}{Averaged Branching Fraction}      \\ \hline
  $D_{s}^{+}\to\phi e^+\nu_e$                        & $(2.37\pm 0.11)\%$& \cite{ref:BESIII2020phiBF,ref:BESIII2018phiBF,ref:f0PDGBF,ref:BABAR2008phiBF} \\ \hline
$D_{s}^{+}\to\eta e^+\nu_e$                          & $(2.32\pm 0.08)\%$& \cite{ref:BESIII2019etaBF,ref:BESIII2016etaBF,ref:f0PDGBF} \\ \hline
$D_{s}^{+}\to\eta'e^+\nu_e$                          & $(0.80\pm 0.07)\%$ &\cite{ref:BESIII2019etaBF,ref:BESIII2016etaBF,ref:f0PDGBF} \\ \hline
$D_{s}^{+}\to K^0 e^+\nu_e$                          & $(0.34\pm 0.04)\%$&\cite{ref:BESIII2019KBF,ref:f0PDGBF} \\ \hline
$D_{s}^{+}\to K^{*}(892)^0 e^+\nu_e$                 & $(0.21\pm 0.03)\%$& \cite{ref:BESIII2019KBF,ref:f0PDGBF} \\\hline
$D_{s}^{+}\to f_{0}(980)e^+\nu_e, \;f_{0}(980)\to \pi\pi$                 & $(0.30\pm 0.05)\%$&\cite{ref:newf0PDGBF}\\ \hline\hline
Sum of Semielectronic Modes                  & $(6.34\pm 0.17)\%$& \\ \hline\hline
$\BF\left(\semilep\right)$                                 & $(6.5 \pm 0.4)\%$& \cite{ref:CLEOMeas}                    \\ \hline\hline
$D_{s}^{+}\to\tau^+\nu_\tau\to e^+\overline{\nu}_\tau\nu_e\nu_\tau$ & $(0.96\pm 0.04)\%$& \cite{ref:BESIII2016tauBF,ref:Belle2013tauBF,ref:BABAR2010tauBF,
ref:CLEO2009taurhoBF,ref:CLEO2009taupiBF,ref:CLEO2009taueBF}     \\\hline
\end{tabular}
}

\raggedright
\caption{Branching fractions for observed $D_{s}^{+}$ semielectronic decays and for $\taucont$. The listed uncertainties are the total uncertainties. The $D_{s}^{+}\to f_{0}(980)e^+\nu_e, \;f_{0}(980)\to \pi\pi$ branching fraction is calculated based on the measurements of  $D_{s}^{+}\to\pi^+\pi^-e^+\nu_e$ from~\cite{ref:newf0PDGBF} and corrected for the corresponding $\pi^0\pi^0$ branching fraction by an isospin correction factor of $\frac{3}{2}$.}
\label{table:BFS}
\end{table*}

\renewcommand{\arraystretch}{1.}

Table~\ref{table:BFS} lists the six exclusive $D_s^+$ semielectronic
modes that have been observed to date and their branching
fractions, as well as the branching fraction of $\taucont$
and the previously measured $D_{s}^{+}\to Xe^{+}\nu_e$ branching fraction
\cite{ref:CLEOMeas,ref:BESIII2020phiBF,ref:BESIII2018phiBF,ref:f0PDGBF,ref:newf0PDGBF,ref:BABAR2008phiBF,ref:BESIII2019etaBF,ref:BESIII2016etaBF,ref:BESIII2019KBF,ref:BESIII2016tauBF,ref:Belle2013tauBF,ref:BABAR2010tauBF,ref:CLEO2009taurhoBF,ref:CLEO2009taupiBF,ref:CLEO2009taueBF}. By
comparing the inclusive semielectronic branching fraction with the sum
of all measured exclusive semielectronic branching fractions, we can
estimate the branching fraction for further unobserved
$D_s^+$ semielectronic decays.  Measurements of the semielectronic
branching fractions for different charmed mesons can be combined to
probe for non-spectator eﬀects in heavy-meson decays
\cite{ref:Bigi,ref:Voloshin}.  The CLEO-c measurement of the ratio of
the $D_s^+$ and $D^0$ semielectronic widths,
$\frac{\Gamma(D_s^+\rightarrow X e^+ \nu_e)}{\Gamma(D^0\rightarrow X
  e^+ \nu_e)} = 0.815 \pm 0.052$, is in agreement with predictions
employing an eﬀective quark model~\cite{ref:Rosner} and shows that
non-spectator eﬀects are present in semielectronic charmed-meson
decays.  It has also been demonstrated with CLEO-c data that the
inclusive semielectronic momentum spectrum can be used to make
sensitive tests for specific non-spectator processes, such as Weak
Annihilation (WA)~\cite{ref:Manohar,ref:Gambino}. Strong understanding of these processes are required for application of heavy-quark-expansion in extracting CKM elements from inclusive semileptonic $B$ meson decays~\cite{ref:Vos}.  Thus, the improved
precision of both the inclusive branching fraction and the momentum
spectrum of $D_s^+$ semielectronic decays reported in this article
have potential to contribute to reducing theoretical uncertainties in
determining CKM parameters with heavy-meson decays.

The remainder of the article is organized in seven sections. The BESIII detector, the analyzed data, and the Monte Carlo (MC) simulation samples are described in Sec.~\ref{sec:detector}. An overview of the measurement technique is presented in Sec.~\ref{sec:technique}. Event-selection requirements based on full reconstruction of hadronic $D_s^-$ decays are discussed in Sec.~\ref{sec:stag}. Semielectronic decay selection requirements and further analysis of candidate signal events are presented in Sec.~\ref{sec:analysis}. The systematic uncertainties of our measurement are evaluated in Sec.~\ref{sec:sys}. We conclude with a summary of our results in Sec.~\ref{sec:conc} and acknowledgements in Sec.~\ref{sec:acknowledgements}.

\section{Detector and data samples}\label{sec:detector}

The BESIII detector records the results of symmetric $e^+e^-$ collisions provided by the BEPCII collider~\cite{ref:BEPCII}. BEPCII produces collisions at center-of-mass energies ($E_{\text{cm}}$) between $2$ and $4.9$ GeV, and BESIII has collected the world's largest data samples near a number of pair-production threshold energies for charmed hadrons. The BESIII detector is composed of the following sub-systems for particle detection and identification: a helium-based multilayer drift chamber (MDC), a plastic scintillator time-of-flight system (TOF), a CsI(Tl) electromagnetic calorimeter (EMC), a 1.0 T superconducting solenoid, and a set of resistive plate chambers for muon identification (MUC). The acceptance of the BESIII detector for charged particles is $93\%$ of the full solid angle. Ionization energy deposits in the MDC are used to reconstruct charged particle tracks and determine particle momenta from the curvature in the magnetic field. The MDC provides a momentum resolution of $0.5\%$ for particles with a momentum of $1$ GeV/$c$. The end cap TOF system was upgraded in 2015 with multi-gap resistive plate chamber technology, providing a time resolution of 60 ps~\cite{ref:TOF}. Measurements of charged particle specific ionization in the MDC ($dE/dx$), the flight time in the TOF, and the energy deposit in the EMC associated with a track are combined to identify particles. A detailed description of the BESIII detector is given in~\cite{ref:BESIII}.

Inclusive MC samples with the equivalent of five times the luminosity of
data and exclusive
MC samples with the equivalent of thirty-five times the luminosity of
data are used to optimize selection criteria, investigate
distributions of signal and background processes, and determine the
efficiencies of our selection criteria. Hadronic simulation samples
are produced using the event generators KKMC~\cite{ref:KKMC} and
EvtGen
\cite{ref:EvtGen,ref:BesEvtGen}. BabayagaNLO~\cite{ref:BabayagaNLO} is
used to produce radiative Bhabha $(e^+e^-\to\gamma e^+e^-)$
samples. Final-state radiation (FSR) of particles is simulated with
PHOTOS~\cite{ref:PHOTOS}. Interaction of the simulated particles with
the detector material is handled by GEANT4 \cite{ref:GEANT4}, which
uses a detailed XML-based description of the detector
geometry~\cite{ref:Bes3Geometry}.

In 2016, BESIII collected $3.19 \text{ fb}^{-1}$ of data at $\Ecm=4.178\GeV$, which provides approximately $6.4 \times 10^6$ $D_{s}^{+}$ mesons primarily through the process $e^+e^-\toD_{s}^{*+}D_{s}^{-}$, with a small contribution from $e^+e^-\toD_{s}^{+}D_{s}^{-}$. We analyze the entirety of this sample, as well as $2.08\text{ fb}^{-1}$ of data collected at center-of-mass energies in the range $\Ecm=4.189-4.219\GeV$ in 2017, and $1.05\text{ fb}^{-1}$ of data collected in the range $\Ecm=4.225-4.230\GeV$ in 2013. These three samples are analyzed separately due to differing detector and running conditions. \iffalse The center-of-mass energies for each data sample are measured with use of the process $e^+e^-\to\mu^+\mu^-$, and the luminosities are determined through measurement of $e^+e^-\to e^+e^-$.\fi  A summary of the data sets with their $E_{\text{cm}}$~\cite{ref:4230EBeam}, integrated luminosities~\cite{ref:4230Lum}, and estimated number of $D_{s}^{+}$ mesons produced is shown in Table~\ref{table:samples}.

\begin{table}[hbt]
\centering{
\begin{tabular}{c|c|c}
$E_{\text{cm}}$ (MeV)& $\int \mathcal L \;dt$ ($\text{pb}^{-1}$)& $N_{D_s} (\times 10^6)$\\ \hline
$ 4178$ & $3189.0 \pm 0.9 \pm 31.9$ & $6.4$ \\  \hline
$4189$   & $526.7 \pm 0.1 \pm 2.2$  &$1.0$\\
$4199$   & $526.0 \pm 0.1 \pm 2.1$  &$1.0$\\
$4209$  & $517.1 \pm 0.1 \pm 1.8$  &$0.9$\\
$4219 $   & $514.6 \pm 0.1 \pm 1.8$  &$0.8$\\ \hline
$4225-4230$ \cite{ref:4230EBeam}  &$1047.3 \pm 0.1 \pm 10.2$ \cite{ref:4230Lum}&$1.3$
\end{tabular}
}

\raggedright
\caption{$E_{\text{cm}}$, integrated luminosities, and estimated number of $D_s$ mesons $(N_{D_s})$ for the analyzed data samples. In each case, the first listed uncertainty is statistical and the second is systematic.}
\label{table:samples}
\end{table}

\section{Measurement Technique}\label{sec:technique}

Our analysis procedure employs the double-tag technique pioneered by the MARKIII collaboration~\cite{ref:MARKIII}. We fully reconstruct hadronic $D_s^-$ mesons (the tag or tag-side), and determine the number of signal events by analyzing the remaining charged tracks unused in the reconstruction of the tags (the recoil or recoil-side). We refer to events where a tag meson is found as single-tag events, and events where a semielectronic decay is identified on the recoil side in addition to the tag meson as double-tag events. Our single-tag selection is described in more detail in Sec.~\ref{sec:stag}.

Since identification and reconstruction of positron tracks is only possible above a certain momentum threshold with the BESIII detector, double-tag candidate tracks are only accepted with $p > 200 \text{ MeV}/c$. The differential inclusive semielectronic branching fraction as a function of the electron momentum $p_e$ ($ \frac{d\mathcal B_\text{SL}}{dp_e}$)  is given as

\begin{equation}
\frac{d\mathcal B_\text{SL}}{dp_e}=\frac{\Delta\left(\frac{ n_\text{DT}}{\epsilon_\text{DT}}\right)/\Delta p_{e}}{n_\text{ST}/\epsilon_\text{ST}}=\frac{\left(\frac{\Delta n_\text{DT}}{\epsilon_\text{sig}}\right)/\Delta p_e}{n_\text{ST}\frac{\epsilon'_\text{ST}}{\epsilon_\text{ST}}}=\frac{\Delta N_\text{DT}/\Delta p_e}{n_\text{ST}b_{\text{tag}}}.
\label{eqn:DiffBRF}
\end{equation}
In this equation $n_\text{ST}$ is the number of the observed single-tag events, $\Delta n_\text{DT}$ is the number of the observed (semielectronic) double-tag events in a particular $p_e$ momentum bin, and $\epsilon_\text{ST}$ and $\epsilon_\text{DT}$ are the single-tag and double-tag reconstruction efficiencies, respectively. We define $\epsilon_\text{DT}=\epsilon'_\text{ST}\times \epsilon_\text{sig}$, where $\epsilon_\text{sig}$ is the momentum-dependent efficiency of reconstructing only the signal side, and $\epsilon'_\text{ST}$ is the tag-side reconstruction efficiency given that a signal event is present on the recoil-side. We define $b_{\text{tag}}\equiv \frac{\epsilon'_\text{ST}}{\epsilon_\text{ST}}$ as the tag bias, which accounts for the difference in the single-tag detection efficiency given that a semielectronic decay is present on the recoil side, and $\Delta N_\text{DT}$ as the true number of double-tag signal events in our single-tag sample for a particular $p_e$ momentum bin.

We begin our signal selection by sorting recoil-side tracks with $p>
200 \text{ MeV}/c$ into momentum bins.  For a given momentum bin, we
further sort these tracks by their charge: tracks with charge opposite
to the tag (right-sign or RS) and tracks with the same charge as the
tag (wrong-sign or WS). The WS sample is used to estimate
charge-symmetric backgrounds in the RS distributions, which include
the true signal. For both charge categories in each momentum bin, we
sort tracks into three mutually exclusive particle identification
hypotheses: electron/positron ($e$ ID), pion ($\pi$ ID), and kaon ($K$
ID). The similar detector response to muons and pions allows us to
treat muons as identical to pions with negligible uncertainty. For
each of these six categories in a given momentum bin, we determine the
number of tracks that originate from a true $D_s^+$ meson by fitting
to the invariant-mass distribution of the tag $D_s^-$ candidate.

We relate the true number of double-tag signal events in our single-tag sample and the observed number of signal-candidate tracks from double-tag events through detector response matrices for charged particle tracking ($A_\text{trk}$) and PID ($A_\text{PID}$):

\begin{equation}\label{eqn:unfolding}
n^b_\text{DT}(p_{i})=\sum_a A_\text{PID}(b|a,p_{i})\sum_{j}A_\text{trk}(a|p_{i},p_{j})N_\text{DT}^a(p_{j}).
\end{equation}
Here, $A_\text{trk}(a|p_{i},p_{j})$ gives the probability of a
particle of type $a$ ($a\in \{e,\pi,K\}$) tracked in momentum bin
$p_{i}$ to have true momentum in bin
$p_{j}$. $A_\text{PID}(b|a,p_{i})$ gives the probability of a particle
of type $a$ passing the PID requirements of particle type $b$ ($b\in
\{e,\pi,K\}$) in track momentum bin $p_{i}$.

We determine the observed double-tag yields $n^b_\text{DT}(p_{i})$ by fitting the invariant-mass distribution of tag-side mesons. For both RS and WS tracks in each momentum bin, we account for $e$ ID efficiencies and particles faking $e$ ID with the PID unfolding method introduced in Eq.~\eqref{eqn:unfolding}:

\begin{equation}
\small
\label{eqn:PIDUnf}
\begin{bmatrix}
n^e_\text{trk}\\
n^\pi_\text{trk}\\
n^K_\text{trk}\\
\end{bmatrix}=A^{-1}_\text{PID}
\begin{bmatrix}
n^e_\text{DT}\\
n^\pi_\text{DT}\\
n^K_\text{DT}\\
\end{bmatrix},
 \;  A_\text{PID}=\begin{bmatrix}
\epsilon_{e} & \text{P}_{\pi\to e} & P_{K\to e} \\
P_{e\to \pi} & \epsilon_{\pi} & P_{K\to \pi} \\
P_{e\to K} & P_{\pi\to K} & \epsilon_{K}
\end{bmatrix}.
\end{equation}
Here, $\epsilon_a$ is the efficiency of our PID selection requirements for the particles of species $a$ and $P_{a\to b}$ is the probability of a particle of species $a$ passing the selection requirements for a particle of species $b$. We apply this PID unfolding to the fitted RS and WS yields for each momentum bin and take the $e$ solution to determine the number of tracked positrons originating from a $D_s^+$ meson. We then take the difference of the number of tracked RS positrons and WS positrons to subtract the contributions of positrons from Dalitz decays of light mesons in the final state of $D_{s}^{+}$ decay such as $\pi^0\to\gamma e^+e^-$ and any other charge-symmetric backgrounds.

This leaves us with the momentum spectrum of tracked positrons from $D_{s}^{+}\to Xe^{+}\nu_e$ events and a smaller contribution from $D_s^+\to\tau^+\nu_\tau\to e^+\nu_e\nu_\tau\overline\nu_\tau$. To account for tracking inefficiency and mismeasurement, we execute the second unfolding from Eq.~\eqref{eqn:unfolding}:

\begin{equation}
\label{eqn:TrkUnf}
N_\text{DT}(p_j)=A^{-1}_\text{trk}n^e_\text{trk}(p_i).
\end{equation}
This provides the true sum of $D_{s}^{+}\to Xe^{+}\nu_e$ and $ \taucont$ positron momentum spectra. We then subtract the contribution of $\taucont$ based on the branching fraction in Table~\ref{table:BFS} to obtain the $D_{s}^{+}\to Xe^{+}\nu_e$ momentum spectrum for $p_{e}>200\text{ MeV}/c$.

To account for positrons with $p_{e}\leq 200 \text{ MeV}/c$, we fit the expected $D_{s}^{+}\to Xe^{+}\nu_e$ momentum spectrum with an assumed spectrum, described in more detail in Sec.~\ref{sec:analysis}, and add the fitted yields in the region $p_{e}\leq 200 \text{ MeV}/c$ to the measured yields with $p_{e}>200\text{ MeV}/c$. This gives us the signal-efficiency-corrected number of signal events
\begin{equation}\label{eq:extrap}
N_\text{DT}=N_{p_e>200 \text{ MeV/} c}+N_{p_e\leq 200 \text{ MeV/} c},
\end{equation} which allows us to calculate the branching fraction similarly to Eq.~\eqref{eqn:DiffBRF}  through

\begin{equation}
\BF\left(\semilep\right)=\frac{N_\text{DT}}{n_\text{ST}b_{\text{tag}}},
\label{eqn:DTagBRF}
\end{equation}
with the MC simulation prediction of $b_{\text{tag}}$ and the determined number of single-tag events $n_\text{ST}$.

\section{Single-Tag Event Selection}\label{sec:stag}

We select single-tag candidates using only the $D_s^-\to K^+K^-\pi^-$
hadronic decay mode because it is unique in having sufficient
statistics and well-known backgrounds. All tag candidate daughters are
required to satisfy the following track-quality requirements: the
track must originate from a region within 1 (10) cm of the $e^+e^-$
interaction point perpendicular (parallel) to the $z$ axis, which is
the symmetry axis of the MDC, and must be tracked with an angle
$\theta$ with respect to the $z$ axis that satisfies $|\cos
\theta|<0.93$. We apply $\pi$/$K$ ID based on $dE/dx$ and TOF
measurements to all tag-candidate daughters to maximize the purity of
the tag sample.  Multiple candidates in a single event are allowed,
both to increase the single-tag selection efficiency and to minimize
the tag bias, $b_{\text{tag}}$. We calculate the recoil mass against the tag
$D_s^-$ candidate,

\begin{equation}
\small
 M_{\text{Rec}} c^2=\left[\left(E_{\text{cm}}-\sqrt{\left|\vec{p}_{D_s}c\right|^2+m_{D_s}^2c^4}\right)^2-\left|\vec{p}_{D_s}c\right|^2\right]^{\frac 12},
 \end{equation}
where $\vec{p}_{D_s}$ is the momentum of the reconstructed tag $D_s^-$ candidate, and $m_{D_s}$ is the known $D_s$ mass~\cite{ref:PDG}.  We require the recoil mass to be consistent with a $D_{s}^{*+}D_{s}^{-}$ event hypothesis by imposing requirements for each data set: $2057 <M_{\text{Rec}}<2177 \text{ MeV}/c^2$ for $\Ecm=4.178\GeV$, $2145 <M_{\text{Rec}}< 2190  \text{ MeV}/c^2$ for $\Ecm=4.189-4.219\GeV$, and $2150 <M_{\text{Rec}}< 2200 \text{ MeV}/c^2$ for $\Ecm=4.225-4.230\GeV$. For about 20\% of cases, our event selection identifies more than one candidate per event after these recoil-mass requirements. Background due to extra candidates is subtracted as part of the fitting procedure described later in this section. The $M_{\text{Rec}}$ distributions for each data set can be seen in Fig.~\ref{fig:mRec}. While differences between the MC simulation and data can be seen in the figure, our measurement is not sensitive to the MC simulation of this distribution.

\begin{figure}[hbt]
\centering{
\begin{tabular}{c}
\includegraphics[width=3.25in]{figures/4180_mRec.pdf}\\\includegraphics[width=3.25in]{figures/XYZ_mRec.pdf}\\\includegraphics[width=3.25in]{figures/4230_mRec.pdf}\\
\end{tabular}
}

\caption{Recoil-mass distribution of $D_s^-\to K^+K^-\pi^-$ candidates from the  $\Ecm=4.178\GeV$ (top), $\Ecm=4.189-4.219\GeV$ (center), and $\Ecm=4.225-4.230\GeV$ (bottom) data sets. The black points are data, and the black and red histograms are the signal and  background distribution predicted by the MC simulation, scaled by integrated luminosity.  The dashed red lines indicate the selection conditions.}
\label{fig:mRec}
\end{figure}

We determine the number of true single-tag candidates by performing an
unbinned fit to the distribution of the invariant mass of the tag
$D_s^-$ candidates, $M_\text{Inv}$. The signal shape is based on MC
simulation and convolved with a Gaussian function whose width and mean
are left free in the fit to account for a possible difference in
resolution and calibration between data and the MC simulation. The
distribution of backgrounds in this variable is modeled using a
second-order Chebyshev polynomial. The fit range is chosen to be
within $\pm40 \text{ MeV}/c^2$ of the known $D_{s}^{+}$ mass for the $\Ecm=4.178\GeV$
distribution corresponding to approximately $5\sigma$ of the simulated
$M_\text{Inv}$ distribution at $\Ecm=4.178\GeV$. In the $\Ecm=4.189-4.219\GeV$ and $\Ecm=4.225-4.230\GeV$
data sets, the fit range is reduced to be within $\pm35 \text{ MeV}/c^2$
of the known $D_s$ mass due to non-polynomial background structures
appearing at the edges of the fit range. The single-tag fits to each
data set are shown in Fig.~\ref{fig:mInv} and the fitted yields are
listed in Table~\ref{table:STYields}. The $\Ecm=4.178\GeV$ data set provides
approximately $56\%$ of our single-tag sample, and the $\Ecm=4.189-4.219\GeV$ and
$\Ecm=4.225-4.230\GeV$ data sets provide $33\%$ and $11\%$, respectively. Within each
data set, $99\%$ of our $D_s^{-}$ candidates come from $D_s^*D_s$
events and $1\%$ come from $D_sD_s$ events. $D_s^*D_s^*$ events can be produced at $\Ecm=4.225-4.230\GeV$, but contribute negligibly to our candidate sample.

The single-tag reconstruction efficiency depends on the topology of the recoil-side decay, which we account for through $b_{\text{tag}}$. We determine the tag bias from MC samples for double-tag events from each of the observed semielectronic modes listed in Table~\ref{table:BFS} and for $\taucont$. The $D_{s}^{+}\to Xe^{+}\nu_e$ tag bias is determined by averaging the tag biases for the six observed modes weighted by the branching fractions in Table~\ref{table:BFS}. The single-tag efficiencies with no specification on the signal-side and the determined tag biases for $D_{s}^{+}\to Xe^{+}\nu_e$ and $\taucont$ for each data set are shown in Table~\ref{table:bias}.
\begin{table}[H]
\centering{
\begin{tabular}{c|c}
Data Set & Fitted Single-Tag Yields\\ \hline
$\Ecm=4.178\GeV$ & $147581\pm 779$\\ \hline
$\Ecm=4.189-4.219\GeV$ & $85845\pm 705$\\ \hline
$\Ecm=4.225-4.230\GeV$ & $29234\pm 435$\\ \hline\hline
Sum &  $262660 \pm 1137$ \\ \hline
\end{tabular}
}

\raggedright
\caption{Fitted single-tag yields from each data set. Only statistical uncertainties are shown.}
\label{table:STYields}
\end{table}

\begin{figure}[h!]
\centering{
\begin{tabular}{ccc}
\includegraphics[width=3.25in]{figures/4180_data_mInv.pdf}\\
\includegraphics[width=3.25in]{figures/XYZ_data_mInv.pdf}\\
\includegraphics[width=3.25in]{figures/4230_data_mInv.pdf}
\end{tabular}
}
\caption{Results of fits to invariant-mass distributions of $D_s^-\to K^+K^-\pi^-$ single-tag candidates from the $\Ecm=4.178\GeV$ data set (top), $\Ecm=4.189-4.219\GeV$ data set (center),  and $\Ecm=4.225-4.230\GeV$ data (bottom). In each plot, the solid blue line is the result of the total fit, the dashed red line is the fitted distribution of non-$D_s^-\to K^+K^-\pi^-$ backgrounds, the dotted black line is the fitted signal distribution, the filled red histogram is the contribution from backgrounds predicted by MC simulation, and the black points are data. The choice of binning is arbitrary and used solely for display. MC simulation distributions are scaled by luminosity.}
\label{fig:mInv}
\end{figure}

\begin{table*}[hbtp]
\centering{
\begin{tabular}{c|c|c|c}
Data Set&$\epsilon_\text{ST}$ &$b_{\text{tag}}$ &$b^\tau_{\rm{tag}}$ \\ \hline
$\Ecm=4.178\GeV$ & $\left(43.10 \pm 0.01\right)\%$ & $1.007 \pm 0.001$ & $1.037 \pm 0.003$ \\ \hline
$\Ecm=4.189-4.219\GeV$ & $\left(42.46 \pm 0.02\right)\%$ & $1.004 \pm 0.002$ & $1.034 \pm 0.004$ \\ \hline
$\Ecm=4.225-4.230\GeV$ & $\left(40.48 \pm 0.03\right)\%$ & $1.005 \pm 0.003$ & $1.044 \pm 0.007$
\end{tabular}
}
\caption{Single-tag efficiencies for all events ($\epsilon_\text{ST}$), and tag bias for $D_{s}^{+}\to Xe^{+}\nu_e$ events ($b_{\text{tag}}$) and for $\taucont$ events ($b^\tau_{\rm{tag}}$) from each data set determined with MC  samples. Listed uncertainties are only statistical.}
\label{table:bias}
\end{table*}

\section{Double-Tag Selection And Analysis}\label{sec:analysis}
  After a single-tag candidate is found, we begin searching for positron candidates among the recoil-side tracks. We sort recoil-side tracks that satisfy track requirements defined in Sec.~\ref{sec:stag} into eighteen 50-MeV$/c$ momentum bins between $200 \text{ MeV}/c$ and $1100 \text{ MeV}/c$. Mutually exclusive PID hypotheses are assigned to each track based on information from the MDC, TOF, and EMC.

Based on our PID assigment, we fill tag-side invariant-mass
distributions in each momentum bin for each of the six categories of
recoil-side tracks defined in Sec.~\ref{sec:technique}. We determine
the number of tracks in each category originating from true $D_{s}^{+}$
events by performing an unbinned fit to the tag-side invariant mass
distribution with the same signal shape used for the single-tag fit,
with the Gaussian parameters fixed to those determined in the
single-tag fit.  We employ a first-order Chebyshev polynomial to model
the distribution of backgrounds. Two examples of the 172 fits we
perform are shown in Fig.~\ref{fig:DatamInvFits}. The full set of fits
is made available as supplemental
material~\cite{ref:supplement}. Measured yields as a function of
momentum bin for each category are shown in Fig.~\ref{fig:yields}.

\begin{figure}[hbt]
\centering{
\begin{tabular}{c}
\includegraphics[width=3.25in]{figures/4180_data_Barrel_450_500_EID.pdf}\\
\includegraphics[width=3.25in]{figures/4180_data_Barrel_450_500_EID_ws.pdf}
\end{tabular}
}

\caption{Example tag invariant-mass fits for double-tag RS (top) and
  WS (bottom) positron candidates with momentum in the range 450-500
  MeV$/c$ from $\Ecm=4.178\GeV$ data. In both plots, the solid blue line is the
  result of the total fit, the dashed red line is the fitted
  distribution of non-$D_s^-\to K^+K^-\pi^-$ backgrounds, the dotted
  black line is the fitted signal distribution, the filled red
  histogram is the contribution from backgrounds predicted by MC
  simulation, and the black points are data. The choice of binning is
  arbitrary and used solely for display.  MC simulation distributions
  are scaled by the number of single-tag events.}
\label{fig:DatamInvFits}
\end{figure}

\begin{figure*}[hbt]
\centering{
\large
\begin{tabular}{cc}
RS $e$ yields& WS $e$ yields\\
\includegraphics[width=3.00in]{figures/EID_yields_data.pdf}&\includegraphics[width=3.00in]{figures/EID_ws_yields_data.pdf}\\
 RS $\pi$ yields& WS $\pi$ yields\\
\includegraphics[width=3.00in]{figures/PiID_yields_data.pdf}&\includegraphics[width=3.00in]{figures/PiID_ws_yields_data.pdf}\\
 RS $K$ yields& WS $K$ yields\\
\includegraphics[width=3.00in]{figures/KID_yields_data.pdf}&\includegraphics[width=3.00in]{figures/KID_ws_yields_data.pdf}
\end{tabular}
}
\caption{Measured yields for each track category as a function of momentum. The results from $\Ecm=4.178\GeV$ are shown as black squares, those from $\Ecm=4.189-4.219\GeV$ as blue triangles, and those from $\Ecm=4.225-4.230\GeV$ as red circles.}
\label{fig:yields}
\end{figure*}

With the yields determined from the fits, we perform the matrix unfolding procedure described in Sec.~\ref{sec:technique} to correct for the inefficiencies of our electron/positron identification and for the misidentification of pions and kaons as electrons in each momentum bin for both RS and WS tracks, as described in Eq.~\eqref{eqn:PIDUnf}. The elements of $A_\text{PID}$ for each data set are determined by applying our PID requirements to MC samples of particles originating from $D_{s}^{+}$ decays. Differing detector conditions of the three data sets introduce deviations among the PID rates from the data sets. The rates that populate the $A_\text{PID}$ matrices are shown in Fig.~\ref{fig:PIDeffs} for each data set. The results of the unfolding procedure for both RS and WS tracks are shown in Fig.~\ref{fig:PIDUnfolding}, as well as the difference of the PID-unfolded right-sign and wrong-sign yields, which gives the PID-unfolded momentum spectra for tracked positrons originating from $D_{s}^{+}\to Xe^{+}\nu_e$ and $\taucont$ events.  The assessment of systematic uncertainties related to the two rates to which we are most sensitive (the electron identification efficiency and the pion-faking-electron rate) are discussed in Sec.~\ref{sec:sys}.

\begin{figure*}[hbt]
\centering{
\large
\begin{tabular}{ccc}
$\epsilon_{e}$&$P_{\pi\to e}$&$P_{K\to e}$\\
\includegraphics[width=2in]{figures/Elec_EID.pdf}&\includegraphics[width=2in]{figures/Pion_EID.pdf}&\includegraphics[width=2in]{figures/Kaon_EID.pdf}\\
$P_{e\to \pi}$&$\epsilon_{\pi}$&$P_{K\to \pi}$\\
\includegraphics[width=2in]{figures/Elec_PiID.pdf}&\includegraphics[width=2in]{figures/Pion_PiID.pdf}&\includegraphics[width=2in]{figures/Kaon_PiID.pdf}\\
$P_{e\to K}$&$P_{\pi\to K}$&$\epsilon_{K}$\\
\includegraphics[width=2in]{figures/Elec_KID.pdf}&\includegraphics[width=2in]{figures/Pion_KID.pdf}&\includegraphics[width=2in]{figures/Kaon_KID.pdf}
\end{tabular}
}

\caption{PID rates as a function of momentum used to populate the $A_\text{PID}$ matrices. Entries from $\Ecm=4.178\GeV$ are shown as black squares, those from $\Ecm=4.189-4.219\GeV$ as blue triangles, and those from $\Ecm=4.225-4.230\GeV$ as red circles.}
\label{fig:PIDeffs}
\end{figure*}

\begin{figure}[!hbt]
\centering{
\begin{tabular}{c}
 \includegraphics[width=3.25in]{figures/4180Data_PIDUnfolding.pdf}\\
 \includegraphics[width=3.25in]{figures/XYZData_PIDUnfolding.pdf}\\
\includegraphics[width=3.25in]{figures/4230Data_PIDUnfolding.pdf}
\end{tabular}
}
\caption{PID unfolding results of positrons from
  $D_{s}^{+}\to Xe^{+}\nu_e$ and $\taucont$ events for the data sets with $\Ecm=4.178\GeV$
  (top), $\Ecm=4.189-4.219\GeV$ (center), and $\Ecm=4.225-4.230\GeV$ (bottom). RS PID-corrected
  yields are shown as green diamonds, WS PID-corrected yields are
  shown as red triangles, their differences are shown as black
  circles, and the predictions from MC samples scaled by the number of
  single-tag events are shown as blue histograms.}
\label{fig:PIDUnfolding}
\end{figure}

After taking the difference of the RS and WS PID-corrected yields, we
apply the tracking unfolding matrix to correct for tracking
reconstruction efficiencies and momentum bin mis-assignment
(Eq.~\eqref{eqn:TrkUnf}). The momentum bin mis-assignment is caused not
only by imperfect detector resolution, but also by FSR from electrons
and positrons, which increases the likelihood for the track momentum
to be less than the momentum produced in the $D_{s}^{+}$ decay. The tracking
unfolding matrices are consistent among the MC samples. \iffalse The
tracking unfolding matrix for the $\Ecm=4.178\GeV$ data set is shown in
Fig.~\ref{fig:TrkUnfMat}.\fi \iffalse The results of the tracking
unfolding procedure are shown in Fig.~\ref{fig:TrkUnfolding}. \fi

We determine the $D_{s}^{+}\to Xe^{+}\nu_e$ positron momentum spectrum for each data sample by subtracting the contribution of positrons from $\taucont$ events from the spectra obtained from the tracking unfolding procedure. We take the $\taucont$ positron momentum spectra from MC samples. We fix the normalization according to

\begin{equation}
N_{DT,\tau}=\mathcal B(\taucont)\times n_\text{ST}\times b_\text{tag,$\tau$},
\end{equation}
where the values of $\mathcal{B}(\taucont)$ and $b_\text{tag,$\tau$}$ are given in Table~\ref{table:BFS}  and Table~\ref{table:bias}, respectively. The subtracted $\tau$ component is shown for each data set in Fig.~\ref{fig:TauSub}.

\begin{figure}[!h]
\centering{
\begin{tabular}{c}
\includegraphics[width=3.25in]{figures/4180Data_TauSub.pdf}\\
\includegraphics[width=3.25in]{figures/XYZData_TauSub.pdf}\\
\includegraphics[width=3.25in]{figures/4230Data_TauSub.pdf}
\end{tabular}
}
\caption{Results of subtracting the contribution of $D_s^+\to\tau^+\nu_\tau\to e^+\nu_e\overline{\nu}_\tau\nu_\tau$ in $\Ecm=4.178\GeV$ data (top), $\Ecm=4.189-4.219\GeV$ data (center), and $\Ecm=4.225-4.230\GeV$ data (bottom). Tracking-corrected yields are shown as blue diamonds, and the predicted contributions from $\taucont$ are shown as red histograms. The differences of the two, the $D_{s}^{+}\to Xe^{+}\nu_e$ yields, are shown as black circles. The dotted black line indicates zero events.}
\label{fig:TauSub}
\end{figure}

By determining the signal-efficiency-corrected number of $D_{s}^{+}\to Xe^{+}\nu_e$ events in each data set with $p_{e}> 200 \text{ MeV}/c$, all effects of the detector response, except for small effects in $b_{\text{tag}}$, have been accounted for. As such, we sum the determined yields from each data set in each momentum bin to produce a combined $D_{s}^{+}\to Xe^{+}\nu_e$ momentum spectrum. A table with the total summed yields as a function of momentum along with their statistical uncertainties is provided in the supplemental material~\cite{ref:supplement}. We determine the number of $D_{s}^{+}\to Xe^{+}\nu_e$ events with $p_{e}\leq 200 \text{ MeV}/c$ by fitting a shape based on MC simulation to the combined yields with $p_{e}>200 \text{ MeV}/c$. The shape is constructed by adding the momentum spectra predicted by MC simulation of the six observed exclusive modes ($\phi e^+\nu_e$, $\eta e^+\nu_e$, $\eta' e^+\nu_e$, $K^0 e^+\nu_e$, $K^*(892)^0 e^+\nu_e$, and $f_0(980) e^+\nu_e$) \iffalse and $D_{s}^{+}\to\tau^+\nu\to e^+\overline{\nu}\nu\nu$\fi in proportion to the branching fractions listed in Table~\ref{table:BFS}. The momentum spectra for the $\eta' e^+\nu$, $K^0 e^+\nu_e$, $K^*(892)^0 e^+\nu_e$, and $f_0(980) e^+\nu_e$ modes are taken from MC samples generated with the ISGW2 model~\cite{ref:ISGW2} for the decay of the $D_{s}^{+}$ meson.
We generate separate MC samples for the two largest modes, $\phi e^+\nu_e$ and $\eta e^+\nu_e$, using simple-pole parameterizations of the respective
decay form factors as functions of the lepton-neutrino system squared
four-momentum $q^2$. The form-factor parameters are taken from
measurements of the BABAR collaboration for $\phi e^+\nu_e$~\cite{ref:BABAR2008phiBF} and the BESIII
collaboration for $\eta e^+\nu_e$~\cite{ref:BESIII2019etaBF}.

The integral of the fitted spectrum from $p_{e}=0 \text{ MeV}/c$ to
$p_{e}=200 \text{ MeV}/c$ is added to the number of
signal-efficiency-corrected $D_{s}^{+}\to Xe^{+}\nu_e$ events with $p_{e}> 200 \text{
  MeV}/c$ as described in Eq.~\ref{eq:extrap}, and the statistical
uncertainty is scaled by the ratio by which the total yield increases due
to this correction. The fit to data with the assumed momentum spectrum can be found in
Fig.~\ref{fig:Data_extrapolation}. As a crosscheck, we also fit to the
data sets separately, and find consistent results. The yields in data
both without and with the correction for the data with $p_{e}<200
\text{ MeV}/c$ are shown in Table~\ref{table:extrapolation}.

\begin{figure}[hbt]
\centering{
\includegraphics[width=3.25in]{figures/extrapolation.pdf}
}
\caption{Momentum spectrum fit used to determine the $D_{s}^{+}$ semielectronic decay yield below $200$ MeV/$c$.  The black points are sums of the final measured $D_{s}^{+}\to Xe^{+}\nu_e$ yields (black circles) for the three data sets from Fig.~\ref{fig:TauSub}, and the solid blue line is result of the fit described in the text.}
\label{fig:Data_extrapolation}
\end{figure}

\begin{table*}[hbt]
\centering{
\begin{tabular}{c|p{4cm}|c}
Sample & \centering{Observed Yields with $p_{e}>200 \text{ MeV}/c$ }&Corrected Yields $(N_\text{DT})$\\ \hline
$\Ecm=4.178\GeV$ &\centering{} $8793 \pm 218$&$9628 \pm 238$ \\ \hline
$\Ecm=4.189-4.219\GeV$ &\centering{} $4695 \pm 174$&$5142 \pm 191$ \\ \hline
$\Ecm=4.225-4.230\GeV$ &\centering{} $1712 \pm 104$&$1874\pm 114$ \\ \hline \hline
Combined & \centering{}$15201 \pm 298$ & $16648\pm 326$
\end{tabular}
}
\caption{Yields in data before and after momentum-extrapolation correction. The ``Combined'' row shows the results of summing together the yields from the three data sets and fitting to the summed distribution. Shown uncertainties are only statistical.}
\label{table:extrapolation}
\end{table*}

Using Eq.~\eqref{eqn:DTagBRF}, we determine $\BF\left(\semilep\right)$ with the momentum-extrapolated number of signal-efficiency-corrected double-tag events from Table~\ref{table:extrapolation}, the observed number of  single-tag events from Table~\ref{table:STYields}, and the tag bias from Table~\ref{table:bias}. The branching fractions determined from each data set independently, their average, and the branching fraction determined from their combination are shown along with the associated statistical uncertainties in Table~\ref{table:SemiLepBF}. The difference of the ``Combined'' result from Table~\ref{table:SemiLepBF}, $\BF\left(\semilep\right)=\left(6.30\pm0.13(\text{stat.})\right)\%$, and the sum of the observed exclusive semielectronic branching fractions from Table~\ref{table:BFS} gives the unobserved $D_{s}^{+}$ semielectronic branching fraction as  $\left(-0.04 \pm 0.21\right)\%$, where the stated uncertainty includes the total uncertainty from the exclusive measurements, but only the statistical uncertainty from the inclusive measurement presented in this article.

\begin{table}[hbt]
\centering{
\begin{tabular}{c|c}
Sample & $\BF\left(\semilep\right)$\\ \hline
$\Ecm=4.178\GeV$ &$\left(6.38\pm0.16\right)\%$ \\ \hline
$\Ecm=4.189-4.219\GeV$ &$\left(5.96\pm0.23\right)\%$ \\ \hline
$\Ecm=4.225-4.230\GeV$ & $\left(6.38\pm0.40\right)\%$ \\ \hline \hline
Combined & $\left(6.30\pm0.13\right)\%$
\end{tabular}
}

\caption{$\BF\left(\semilep\right)$ determined from data. The ``Combined'' row shows the results of summing together the yields from the three data sets and fitting to the summed distribution. Shown uncertainties are only statistical.}
\label{table:SemiLepBF}
\end{table}

\section{SYSTEMATIC UNCERTAINTY}\label{sec:sys}

Our methods to determine the relative systematic uncertainty on our measured $\BF\left(\semilep\right)$ are described below.

\subsection{MC Simulation Statistics and Matrix-Inversion Stability}
We probe the effects of finite MC sample statistics and the stability of the matrix-inversion algorithm by creating a toy ensemble of variations for each efficiency matrix used in the analysis (54 PID matrices and three tracking matrices). For each matrix, we create this ensemble by sampling each entry within the MC sample's statistical uncertainty. For each matrix in each ensemble, we perform the matrix inversion and reperform the analysis with the new inverted matrix. All variations produce a negligible change in our final result, which indicates a negligible systematic uncertainty from the statistical uncertainty of the MC samples as well as the stability of the algorithm for inverting our efficiency matrices.

\subsection{Tracking}\label{sec:TrkSys}

Simulation of our tracking efficiency is studied with a control sample of radiative Bhabha events. Tracking efficiencies as a function of momentum are measured in each data set, as well as in MC samples produced with the BabayagaNLO package~\cite{ref:BabayagaNLO}. The ratios of the measured efficiencies in data and MC samples are weighted by the predicted momentum distribution from signal MC simulation and the number of single-tag events in each data set to determine the systematic uncertainty. This results in a relative systematic uncertainty of $0.7\%$.

In addition, we investigate the systematic uncertainty in the individual tracking efficiency matrix entries. As we assign a systematic uncertainty for the total tracking efficiency, we probe the uncertainty in the individual entries by keeping the sum of a row of the matrix constant while varying the individual entries. The specific variation is as follows:
\begin{itemize}
\item $\left(A_\text{trk}\right)_{i,i-1}\Rightarrow 1.5\left(A_\text{trk}\right)_{i,i-1}$
\item $\left(A_\text{trk}\right)_{i,i+1}\Rightarrow 1.25\left(A_\text{trk}\right)_{i,i+1}$
\item $\left(A_\text{trk}\right)_{i,i}$ decreases to keep the sum of the column constant.
\end{itemize}
This variation is chosen as a conservative estimation of the uncertainties from FSR and detector resolution. We see negligible change when we perform such a variation, so we  only assign the previously stated uncertainty for tracking.
\subsection{PID}\label{sec:PIDSys}

Similar to our procedure in assessing the systematic uncertainty in
our tracking efficiencies, we measure $e$ ID efficiencies as a
function of momentum and track angle in radiative Bhabha control
samples for each data set.

We also probe the accuracy of the pion-faking-electron rates from MC simulation via a control sample of pions collected in each data sample through the decay chains $D^{*+}\to\pi^+ D^0,$  $D^0\to K^-\pi^+,K^-\pi^+\pi^+\pi^-$. To determine the total uncertainty from PID rates, we simultaneously vary the $e$ ID efficiencies and the pion-faking-electron rates using the central values of the measured data-to-MC efficiency ratios and reperform our analysis. This yields a $0.8\%$ change in the final branching fraction, which we assign as the relative systematic uncertainty due to PID.

As our sensitivity to kaon-faking-electron rates is small due to the relatively few number of kaons, the systematic uncertainty in kaon-faking-electron rates is neglected.

\subsection{Tag Bias}
We follow the procedure laid out in~\cite{ref:TagBiasRef}, which assigns a fraction of $1-b_{\text{tag}}$ as the systematic uncertainty based on the particles in the final state of the single-tag $D_s^-$ decay. The specific guidelines for variation of detector-response parameters are as follows: $1.0\%$ per kaon for tracking, $0.5\%$ per pion for tracking, and $0.5\%$ per kaon or pion for PID. For $D_s^-\to K^+K^-\pi^-$, with two kaons and one pion, the quadrature sum is $2.9\%$. With $b_{\text{tag}}$ from Table~\ref{table:bias} (including the contribution from $\taucont$), taking $2.9\%$ of $1-b_{\text{tag}}$ yields a $0.03\%$ relative systematic uncertainty. We additionally propagate the uncertainties in the branching fractions (Table~\ref{table:BFS}) through the calculation of the weighted-average $b_{\text{tag}}$. This yields a $0.07\%$ relative systematic uncertainty in the bias. We add these in quadrature and assign the relative systematic uncertainty due to tag bias as $0.1\%$.

\subsection{Number of Single Tags}
We investigate the systematic uncertainty in the invariant-mass fitting procedure used to determine the number of single tags by varying the choice of background distribution from the nominal second-order Chebyshev polynomial. We use both first-order and third-order Chebyshev polynomials as variations in fitting to each data set. Using the first-order Chebyshev polynomial gives a larger difference in the yields in all cases, while not significantly degrading the quality of the fit. We take an average of the changes for each data set weighted by the single-tag yields to determine the systematic uncertainty, which results in a 0.6\%  relative systematic uncertainty in the number of single tags.

\subsection{Background Shapes}
To assess the uncertainty due to our chosen background shapes in our signal-side fits of the tag invariant mass, we use background distributions based on MC simulation instead of the nominal first-order Chebyshev polynomial to model backgrounds in each of the signal-side fits.  We then reperform the analysis with the yields determined from these alternative fits. The relative difference in $N_\text{DT}$ is 0.4\%, which we assign as the relative systematic uncertainty due to this source.

\subsection{$\mathcal{B}(D_{s}^{+}\to\tau^+\nu_\tau)$}
As systematic uncertainty in the kinematic distributions of $\taucont$ events contribute negligible uncertainty, we account for uncertainty in the contribution of $\taucont$ events by propagating the absolute uncertainty for $\mathcal{B}(\taucont)$ from Table~\ref{table:BFS}, $0.04\%$, which yields a relative systematic uncertainty of $0.6\%$ on our measurement of $\BF\left(\semilep\right)$.

\subsection{Spectrum Extrapolation}\label{sec:SpectrumSys}

We assess the uncertainty due to the momentum spectrum extrapolation
by generating an ensemble of alternative momentum spectra and fitting
these to the data. Each spectrum in this ensemble is created by
Gaussian sampling the branching fractions of the six observed
exclusive semielectronic modes and adding spectra for unobserved decay
modes. We then add these spectra in proportion based on the sampled
branching fractions. \iffalse The individual spectra are shown in
Fig.~\ref{fig:spectra}. \fi

We consider effects from the combination of three unobserved decay modes: $D_{s}^{+}\to h_1\left(1415\right)e^+\nu_e$, $D_{s}^{+}\to f_1\left(1510\right)e^+\nu_e$, and $D_{s}^{+}\to\gamma e^+\nu_e$. MC samples for  $D_{s}^{+}\to h_1\left(1415\right)e^+\nu_e$ and $D_{s}^{+}\to f_1\left(1510\right)e^+\nu_e$ are generated based on the ISGW2 model's predictions~\cite{ref:ISGW2}. We generate MC samples for $D_{s}^{+}\to\gamma e^+\nu_e$ based on the model of Yang and Yang~\cite{ref:GammaENuModel}. We determine the normalization of the $D_{s}^{+}\to h_1\left(1415\right)e^+\nu_e$ and $D_{s}^{+}\to f_1\left(1510\right)e^+\nu_e$ spectra by fixing the relative branching fraction of these decays to the ISGW2 predictions and fitting them in addition to our nominal momentum spectrum. The $D_{s}^{+}\to\gamma e^+\nu_e$ spectrum is fixed to its measured $90\%$ confidence level upper limit, $1.3\times 10^{-4}$~\cite{ref:GammaENu}. From this fit, we determine  an upper limit at the $90\%$ confidence-level for the sum of the branching fractions of $D_{s}^{+}\to h_1\left(1415\right)e^+\nu_e$ and $D_{s}^{+}\to f_1\left(1510\right)e^+\nu_e$. In our toy ensemble, their summed spectra are fixed based on this upper limit. We performed this same procedure excluding the $D_{s}^{+}\to\gamma e^+\nu_e$ spectrum, but found the combination of modes to produce the largest variation.

The resulting systematic uncertainty is determined by filling a distribution of the relative change in $\BF\left(\semilep\right)$ between the alternative momentum spectra and our nominal spectrum. The linear sum of the mean ($0.33\%$) and RMS ($0.34\%$) of this distribution is taken as the uncertainty, which gives a $0.7\%$ relative systematic uncertainty.

We also probe uncertainty in the models we employ by using the ISGW2-predicted momentum spectra for the $D_{s}^{+}\to\phi e^+\nu_e$ and $D_{s}^{+}\to\eta e^+\nu_e$ modes. Using these alternative spectra instead of the nominal spectra gives a less than $0.1\%$ difference in the measured branching fraction, so we conclude that any uncertainty due to model dependence in our analysis is negligible.

\subsection{Summary of Systematic Uncertainties}
The assigned relative systematic uncertainties for all sources are listed in Table~\ref{table:systematics}. Our systematic uncertainty is not dominated by any single source, but the largest contributions come from the momentum spectrum extrapolation and imperfect simulation of PID and tracking efficiencies. The total relative systematic uncertainty is obtained from the quadrature sum of the assigned relative uncertainties. This gives a total relative systematic uncertainty of  $1.6\%$.

\begin{table}[hbt]
\centering{
\begin{tabular}{c|c}
Source&Relative Uncertainty\\ \hline
Tracking & 0.7\%\\ \hline
PID & 0.8\%\\ \hline
Spectrum Extrapolation & 0.7\%\\ \hline
Background Shapes & 0.4\%\\ \hline
Number of Tags & 0.6\%\\ \hline
Tag Bias & 0.1\%\\ \hline
$\mathcal{B}(D_{s}^{+}\to\tau^+\nu_\tau)$ & 0.6\%\\ \hline\hline
Total & 1.6\%
\end{tabular}
}
\caption{Systematic uncertainties in the measurement of $\BF\left(\semilep\right)$.}
\label{table:systematics}
\end{table}

\section{Summary and Discussion}\label{sec:conc}

Using data collected by the BESIII detector in the center-of-mass-energy range of $E_{\text{cm}}=4.178-4.230\text{ GeV}$, we measure the inclusive semielectronic branching fraction of the $D_{s}^{+}$ meson to be
\[
\BF\left(\semilep\right)=\left(6.30\pm0.13\;(\text{stat.})\pm 0.10\;(\text{syst.})\right)\%.
\]
We also measure the lab-frame momentum spectrum of the positrons produced in this decay, which can be seen in Fig.~\ref{fig:Data_extrapolation}.

Our result is consistent with the measurement from the CLEO-c experiment\cite{ref:CLEOMeas},
\[
\BF\left(\semilep\right) =(6.52\pm 0.39\;(\text{stat.}) \pm 0.15\;(\text{syst.}))\%,
\]
with a factor of $3$ reduction of the statistical uncertainty and a factor of $1.5$ reduction of the systematic uncertainty. The total precision of our measurement is $2.6\%$,  which corresponds to approximately a $2.5$ times improvement in the total precision compared to the measurement from CLEO-c.

By taking the difference between our measurement of $\BF\left(\semilep\right)$ and the sum of the best available measurements for the exclusive semielectronic modes in Table~\ref{table:BFS}, we calculate the unobserved semielectronic branching fraction to be
\[
\begin{aligned}
\BF\left(\semilep\right)-\sum_i \mathcal B\left(D_{s}^{+}\to X_i e^+\nu_e\right)=\\\left(-0.04\pm 0.13\;(\text{stat.})\pm 0.20\;(\text{syst.})\right)\%,
\end{aligned}
\]
where the systematic uncertainty includes the total uncertainty on the measured exclusive branching fractions. Our measurement provides no evidence for the existence of unobserved $D_{s}^{+}$ semileptonic decay modes and constrains the branching fractions of all unobserved decay modes. In addition, the measured momentum spectrum can be used to further constrain the decay rates of modes with characteristic momentum spectra. The spectrum is included in tabular form in the supplemental material~\cite{ref:supplement}.

With our updated measurement of the $D_{s}^{+}$ semielectronic branching fraction, the CLEO-c measurement of the $D^{0}$ semielectronic branching fraction~\cite{ref:CLEOMeas}, and the 2020 PDG values for the $D_{s}^{+}$ and $D^{0}$ lifetimes~\cite{ref:PDG}, we find
\[
\frac{\Gamma(D_{s}^{+}\rightarrow Xe^+\nu_e)}{\Gamma(D^0\rightarrow Xe^+\nu_e)}=0.790\pm 0.016\;(\text{stat.})\pm0.020\;(\text{syst.}),
\]
where the systematic uncertainty includes the total uncertainty from $\mathcal B \left(D^{0}\to X e^+\nu_e\right)$. This result is in agreement with the prediction of $\frac{\Gamma(D_{s}^{+} \to Xe^+\nu)}{\Gamma(D^{0} \to Xe^+\nu) }=0.813$ from~\cite{ref:Rosner}, supporting the conclusion that the difference in the semileptonic decay widths of $D_{s}^{+}$ and $D^{0}$ mesons can be accounted for within the Standard Model by non-spectator interactions. Further theoretical analysis of our measured spectrum, similar to those of~\cite{ref:Manohar} and \cite{ref:Gambino}, can constrain specific processes like WA of the constituent $c$ and $\overline s$ quarks of the $D_{s}^{+}$, with potential extensions to determinations of $|V_{ub}|$ in semileptonic $B$ decays~\cite{ref:Vos}.

\vspace{3ex}
\section{Acknowledgments}\label{sec:acknowledgements}

The BESIII collaboration thanks the staff of BEPCII and the IHEP computing center for their strong support. This work is supported in part by National Key Research and Development Program of China under Contracts Nos. 2020YFA0406400, 2020YFA0406300; National Natural Science Foundation of China (NSFC) under Contracts Nos. 11625523, 11635010, 11735014, 11822506, 11835012, 11935015, 11935016, 11935018, 11961141012; the Chinese Academy of Sciences (CAS) Large-Scale Scientific Facility Program; Joint Large-Scale Scientific Facility Funds of the NSFC and CAS under Contracts Nos. U1732263, U1832207; CAS Key Research Program of Frontier Sciences under Contract No. QYZDJ-SSW-SLH040; 100 Talents Program of CAS; INPAC and Shanghai Key Laboratory for Particle Physics and Cosmology; ERC under Contract No. 758462; European Union Horizon 2020 research and innovation programme under Contract No. Marie Sklodowska-Curie grant agreement No 894790; German Research Foundation DFG under Contracts Nos. 443159800, Collaborative Research Center CRC 1044, FOR 2359, FOR 2359, GRK 214; Istituto Nazionale di Fisica Nucleare, Italy; Ministry of Development of Turkey under Contract No. DPT2006K-120470; National Science and Technology fund; Olle Engkvist Foundation under Contract No. 200-0605; STFC (United Kingdom); The Knut and Alice Wallenberg Foundation (Sweden) under Contract No. 2016.0157; The Royal Society, UK under Contracts Nos. DH140054, DH160214; The Swedish Research Council; U. S. Department of Energy under Contracts Nos. DE-FG02-05ER41374, DE-SC-0012069.


\section{Introduction}
This very welcome paper by Jacob, Gong, Edlefsen and Dempster \cite{jacob2019gibbs} is a fascinating contribution to an important subject. Persistent questions about the lower and upper probabilities include:
\begin{enumerate}
	\item What does it do in simple problems?
	\item What does it do in real problems?
	\item What do the upper and lower probabilities mean?
	\item How can the needed computations be carried out?
\end{enumerate}

The paper offers answers to these questions by proposing a Gibbs sampler to perform statistical inference for categorical distributions using the Dempster-Shafer approach.  To be precise, let $\boldsymbol{x} = (x_i)_{i=1}^N$ be the observations, each comes from one of the $K$ categories. The model assumes that there exists a $\theta = (\theta_1, \cdots, \theta_K)$ in the $K-$simplex $\Delta :=\{(\theta_1, \theta_2, \cdots, \theta_K): \theta_i \geq 0 \text{~for every~} i, \sum_{i=1}^K \theta_i = 1\}$ such that $\mathbb P(x_i = k) = \theta_k$ for every $i, k$. Moreover, it is  assumed that the observations $\boldsymbol{x}$ are generated from the following procedure.
\begin{enumerate}
	\item For each $i \in \{1, \cdots, N\}$, sample $u_i$ uniformly from the simplex $\Delta$,
	\item Since $\theta$ partitions $\Delta$ into $K$ disjoint pieces, we set $x_i = k$ if $u_i \in \Delta_k(\theta)$, where $\Delta_k(\theta)$ is the ``subsimplex''
	which have the same vertices as $\Delta$ except for the $k$-th vertex replaced by $\theta$.
\end{enumerate}

 Since this machinery generates $x_i$ through $u_i$, it is natural to ask what are the feasible points $\boldsymbol{u} \in  \{(u_1, u_2,\cdots u_n )\in \Delta^N\}$ that can generate the observed $\boldsymbol{x}$. For a fixed $\theta$, the set of feasible points $\mathcal R_{\boldsymbol{x}}(\theta)$ are easy to describe, as it is  the product of $N$ subsimplexes:
$$
\mathcal R_{\boldsymbol{x}}(\theta) = \prod_{i = 1}^N \Delta_{x_i}(\theta),
$$
and the whole feasible set is then
$$\mathcal R_{\boldsymbol{x}} := \cup_{\theta\in \Delta}\mathcal R_{\boldsymbol{x}}(\theta).$$

From $u\in \mathcal R_{\boldsymbol{x}}$ one forms $\mathcal F_u := \{\theta\in \Delta: u_n \in \Delta_{x_n}(\theta) \text{~ for every~} n\}$ and assigns a lower probability to sets $\Sigma$ in the parameter space by the chance that $\mathcal F(u)$ is in $\Sigma$ (this last by Monte Carlo).

The paper's main contribution is a Gibbs sampler that samples from the uniform distribution of $R_{\boldsymbol{x}}$. The problem is challenging due to the complicated nature of  $R_{\boldsymbol{x}}$. In fact,  uniformly sampling from $\mathcal R_{\boldsymbol{x}}(\theta)$ is easy (using Algorithm 1 of \cite{jacob2019gibbs}), but uniform sampling from the union is much more difficult. The key observations in the paper are Proposition 3.1 and Proposition 3.2. They directly characterize the conditional distribution of $\boldsymbol{u}$ and make it possible for implementing the Gibbs sampler. We find the algorithm elegant and insightful.

\section{The math problem}

The new algorithm is interesting as a mathematics problem. The authors have translated it into a clever probability problem in the case of two categories. This gives the Markov chain:
\begin{align}\label{eqn:evolve}
Z^{(t)} = B_1^{(t)} (1 - B_2^{(t)}) Z^{(t-1)} + B_2^{t},
\end{align}
where $B_1^{t} \sim \operatorname{Beta}(N_1, 1)$ and $B_2^{t} \sim \operatorname{Beta}(1, N_2)$ are independent Beta random variables and $N_1, N_2$ are two fixed positive integers.

The chain falls into the well-studied area of iterated random functions  \cite{chamayou1991explicit}\cite{diaconis1999iterated}. It has been previously studied in \cite{letac2002donkey} who offer higher dimensional versions which might be relevant to the present paper when $k\geq 3$.

In the $k= 2$ case, let $P$ be the Markov transition kernel of the chain  and $\pi$ be its stationary distribution. The authors use coupling techniques to derive the following convergence  bound in terms of the Wasserstein-1 distance
\begin{align}\label{eqn: was upper bound}
\mathcal W_1 (P^t(z, \cdot), \pi) \leq \bigg(\frac{N_1}{N_1 + 1} \cdot \frac{N_2}{N_2 + 1}\bigg)^t \cdot W_1 (\delta_{z}, \pi)
\end{align}
for every initialization $z$, where $\delta_z$ is the delta-mass at  $z$. As a small contribution to the conversation, we derive a lower bound of the convergence speed. Our results suggest $\frac{N_1}{N_1 + 1} \cdot \frac{N_2}{N_2 + 1}$ is the exact convergence rate for the Markov chain. To start with, we prove the stationary distribution of the chain is another Beta distribution (this result can also be derived directly from the original formulation of the problem, see Appendix C of \cite{jacob2019gibbs}). The following lemma is well known and helpful:
\begin{lem}\label{lem: beta distribution}
	Let $X\sim \operatorname{Beta}(a,b)$, $Y\sim \operatorname{Beta}(a+b, c)$ be independent random variables with $a,b,c >0$. Then $XY\sim \operatorname{Beta}(a, b+c)$.
\end{lem}

\begin{proof}
	We show the $k$-th moment of $XY$ equals the $k$-th moment of a $\operatorname{Beta}(a, b+c)$ random variable for every $k$.  Since both $XY$ and beta distributions are bounded by $[0,1]$, standard results in probability show that the distribution is characterized by all its moments.

	The $k$-th moment for a $\operatorname{Beta}(a_1, a_2)$ random variable is $\frac{\Gamma(a_1 + k)}{\Gamma(a_1)} \frac{\Gamma(a_1 + a_2)}{\Gamma(a_1 + a_2 + k)}$, therefore,
	\begin{align*}
	\mathbb E(XY)^k &= \frac{\Gamma(a + k)}{\Gamma(a)} \frac{\Gamma(a + b)}{\Gamma(a + b + k)} \times \frac{\Gamma(a + b+ k)}{\Gamma(a + b)} \frac{\Gamma(a + b + c)}{\Gamma(a + b + c + k)} \\
	& =  \frac{\Gamma(a + k)}{\Gamma(a)}  \frac{\Gamma(a + b + c)}{\Gamma(a + b + c + k)},
	\end{align*}
	which is the same as the $k$-th moment of $\operatorname{Beta}(a+b,c)$ distribution, as desired.
\end{proof}

Now we are ready to show $\operatorname{Beta}(N_1+1, N_2)$ is stable under the transformation $Z \rightarrow B_1 (1 - B_2) Z + B_2$, and is thus the stationary distribution of the chain $\{Z^{(0)}, Z^{(1)}, \cdots, Z^{(n)}\}$.
\begin{prop}
	Let $Z\sim \operatorname{Beta}(N_1 + 1, N_2), B_1 \sim \operatorname{Beta}(N_1, 1), B_2 \sim \operatorname{Beta}(1, N_2)$ be independent random variables, then $B_1 (1 - B_2) Z + B_2 \sim \operatorname{Beta}(N_1 + 1, N_2)$.
\end{prop}
\begin{proof}
	Let $Y := B_1 (1 - B_2) Z + B_2$, then  $1 - Y = (1-B_2) (1 - B_1Z)$. Lemma \ref{lem: beta distribution} shows $B_1 Z$ is a $\operatorname{Beta}(N_1, N_2 + 1)$ random variable, which in turn shows $1 - B_1 Z \sim \operatorname{Beta}(N_2+1, N_1)$. Since $1 - B_2 \sim \operatorname{Beta}(N_2, 1)$, apply Lemma \ref{lem: beta distribution} again yields $1- Y \sim \operatorname{Beta}(N_2, N_1 + 1)$. Therefore $Y$ is a $\operatorname{Beta}(N_1 + 1, N_2)$ random variable.
\end{proof}

Now we turn to lower bound the convergence rate of the Markov chain. Recall that for any two probability measure $\mu, \nu$ with bounded support , the Kantorovich- Rubinstein dual theorem  shows:
\[
\mathcal W_1(\mu, \nu)  = \sup_{\{f: \text{Lip}(f)\leq 1\}} \int f (d\mu - d\nu)
\]
Choosing $f(x)  = \pm x$ immediately implies $\mathcal W_1(\mu, \nu) \geq \lvert m_1(\mu) - m_1(\nu)\rvert$, where $m_1$ stands for the first moment. Let $\nu = \pi$ and $Z^{(t)}(z)$ be the Markov chain at time $t$ which starts at $z$. It is clear that $m_1(\nu) = \frac{N_1 + 1}{N_1 + N_2 + 1}$, and the evolution equation \ref{eqn:evolve} yields
\begin{align}\label{eqn:recursion}
\mathbb E(Z^t(z)) = \frac{N_1 N_2}{(N_1 + 1)(N_2 + 1)}  \mathbb E(Z^{t-1}(z)) + \frac{1}{N_2 + 1}
\end{align}
Using the initial condition $\mathbb E(Z^0(z)) = z $, recursion \ref{eqn:recursion} can be solved as:
\[
\mathbb E(Z^t(z))  = \bigg( \frac{N_1 N_2}{(N_1 + 1)(N_2 + 1)}\bigg)^t \bigg(z - \frac{N_1 + 1}{N_1 + N_2 + 1}\bigg) + \frac{N_1 + 1}{N_1 + N_2 + 1}.
\]
Thus we immediately have
\begin{align}\label{eqn: was lower bound}
\mathcal W_1(P^t(z,\cdot),\pi) \geq \bigg\lvert z - \frac{N_1 + 1}{N_1 + N_2 + 1}\bigg \rvert \bigg(\frac{N_1 + 1}{N_1 + N_2 + 1}\bigg)^t.
\end{align}
Therefore, unless the chain is not initialized at exactly $\frac{N_1 + 1}{N_1 + N_2 + 1}$, the exact convergence rate of the Markov chain equals $\frac{N_1 + 1}{N_1 + N_2 + 1}$.
Combining \ref{eqn: was upper bound} and \ref{eqn: was lower bound}, we have the following:
\begin{prop}
	\begin{align}\label{eqn: was convergence rate}
  \lvert\mathbb E Z - z \rvert	\leq \mathcal W_1 (P^t(z, \cdot), \pi)\bigg/ \bigg(\frac{N_1}{N_1 + 1} \cdot \frac{N_2}{N_2 + 1}\bigg)^t \leq \mathbb E\lvert Z - z \rvert,
	\end{align}
	where $Z\sim \operatorname{Beta}(N_1 + 1, N_2)$.
\end{prop}
When $z = 0$ or $1$, Formula \ref{eqn: was convergence rate} shows the convergence speed can be calculately exactly as the lower bound matches the upper bound perfectly. Moreover, the function $u: z \rightarrow \mathbb E|Z - z| = \int_{0}^z \mathbb P (Z \leq t) dt + \int_{z}^1 \mathbb P (Z \geq t) dt$ has derivative $u'(z) = \mathbb P(Z \leq z) -\mathbb P(Z \geq z)$ almost everywhere. Therefore $u$ is first decreasing and then increasing on the unit interval, thus $\max_{z\in[0,1]} u(z) = \max\{u(0), u(1)\}$. This gives us the exact convergence speed for the chain under the worst-case sceneario:
\begin{prop}
\[\sup_{z\in[0,1]}\mathcal W_1 (P^t(z, \cdot), \pi) = \frac{\max\{N_1+1, N_2\}}{N_1 + N_2 + 1} \bigg(\frac{N_1 N_2 }{(N_1 +1)(N_2 + 1)}\bigg)^t.\]
\end{prop}
\begin{proof}
	Since $\mathbb E|Z - z|$ is maximized at either $z = 1$ or $z= 0$, the RHS of Formula $\ref{eqn: was convergence rate}$ is upper bounded by  $\frac{\max\{N_1+1, N_2\}}{N_1 + N_2 + 1}$. Meanwhile, no matter $z = 1$ or $z= 0$, we know $|\mathbb E Z- z| = \mathbb E|Z - z|$, and thus we conclude
	\[\sup_{z\in[0,1]}\mathcal W_1 (P^t(z, \cdot), \pi) = \frac{\max\{N_1+1, N_2\}}{N_1 + N_2 + 1} \bigg(\frac{N_1 N_2 }{(N_1 +1)(N_2 + 1)}\bigg)^t.\]
\end{proof}

\begin{q}
Can the authors see if their general algorithm can be translated into a vector-valued version of \ref{eqn:evolve} for $k$ greater than or equal to $3$?
\end{q}

\section{Hypothesis Testing Problems}
The introduction to the paper emphasizes problems where the sample size $N$ is small compared to the number of categories. We have encountered such problems in our recent work and we find straightforward Bayesian and frequentist solutions.

\begin{eg}
 In studying the popular `wash shuffle' where a deck of $52$ cards is `smushed around the table for $t$ seconds' (say $t = 60$) one wants to test if the cards are well mixed. Here the data consists of $100$ permutations of $X_1, X_2, \cdots X_{100}$ in $S_{52}$. The number of categories $k = \lvert S_{52}\rvert \approx 8.6 \times 10^{68}$ is huge. One has some partial prior information. For example, if the cards are not mixed, it may be because several cards original together are still together. Or the original top (or bottom) card may still be close to the top (or bottom). In our analysis, we found statistics $T_1, T_2, \cdots T_l$ (e.g., $T_1 =$ the number of adjacent pairs, $T_2 =$ position of the original top card, $T_3 = $ the distance from the starting distribution in some natural metric, $\dots$).

 The distribution of these features can be obtained under the null distribution. This allows standard frequentist tests. In our work \cite{diaconis2018bayesian}, we provided a Bayesian solution to this problem. We considered an exponential family through the statistics. For $\sigma \in S_{52}$:
 \[
 \mathbb P_\theta(\sigma) = Z^{-1}(\sigma) \exp(\sum_{i=1}^l \theta_i T_i(\sigma)),
 \]
 where $Z(\sigma) = \sum\limits_{\sigma\in S_{52}} \exp(\sum_{i=1}^l \theta_i T_i(\sigma))$ is the normalization constant.  Thus $\theta = 0$ is the uniform distribution. The exchange algorithm \cite{murray2006mcmc} allowed us to compute the posterior in a reasoning fashion.
\end{eg}

\begin{q}
	Do the authors think there will be a time for their algorithm can be run for such $N$ and $k$?
\end{q}

\section{Remark}
As clearly explained in the paper, the uniform distribution on the $k$-simplex underlies the basic procedure. It is well to be reminded that it is hard to understand the properties when $k$ is large. The following two examples are drawn from \cite{diaconis2002bayesian}.
\begin{eg}[Bayesian Birthday Problem]
	Consider $N$ balls dropped into $k$ boxes, with probability Multinomial$(N,\boldsymbol{\theta})$. What is the chance that all the balls are dropped into distinct boxes?

	If $k = 365$, this is the classical birthday problem. The classical frequentist answer takes $\boldsymbol{\theta} = (\frac1{365}, \cdots, \frac 1{365})$ and shows that the chance is approximately $0.5$ when $N = 23$.  A `flat prior Bayesian' puts a uniform prior on the $k$-simplex. Then the chance of all balls in distinct boxes  is approximately $0.5$ when $N = 16$. Here, being a Bayesian does not change things much.
\end{eg}

\begin{eg}[Bayesian Coupon Collectors Problem]
With notation as above, consider the question: How large does $N$ have to be so that the chance that all boxes are covered is close to $0.5$?

The frequentist solution assumes $\boldsymbol{\theta} = (\frac1{365}, \cdots, \frac 1{365})$ and then shows $N$ should be approximately $2287$. Using a uniform prior for $\boldsymbol{\theta}$, one finds $N$ has to be $191844$. Here, the uniform prior makes a huge difference.
\end{eg}

\begin{q}
	Is the uniform distribution on the simplex a crucial part of the procedure or could this be varied?
\end{q}

\section{A Bit of History}
When P.D. was a beginning graduate student at Harvard (1971), Art Dempster called him in to suggest a possible thesis topic: `Find ways to do the computations required for upper and lower probabilities'. This did not work out at the time but triggered a lifetime's interest. It is inspiring to have tracked his efforts over a $50$ year period. The present paper is important progress. We look forward to progress on the problems mentioned in our introduction.

\bibliographystyle{alpha}