# Carbon Abundance of Globular Cluster M22 (NGC 6656) and the Surface Carbon
Depletion Rates of the Milky Way Globular Clusters

Jae-Woo Lee Department of Physics and Astronomy, Sejong University, 209
Neungdong-ro, Gwangjin-Gu, Seoul 05006, Republic of Korea,
<EMAIL_ADDRESS><EMAIL_ADDRESS>
###### Abstract

It is well known that metal-poor red giant branch (RGB) stars show variations
in some elemental abundances, including carbon, due to the internal mixing
accompanied by their own in situ CN cycle in the hydrogen burning shell. With
our new photometric carbon abundance measurements of RGB stars in M22 and
other globular clusters (GCs) in our previous studies, M5, M3, and M92, we
derive the carbon depletion rates against the $V$ magnitude,
$d\mathrm{[C/Fe]}/M_{V}$, for individual populations in each GC. We find the
metallicity dependence of the carbon depletion rates, $d\mathrm{[C/Fe]}/M_{V}$
$\propto$ $-$0.25[Fe/H]. Our results also suggest that the carbon depletion
rates of the second generation (SG) of stars are larger than those of the
first generation (FG) of stars in our sample GCs, most likely due to different
internal temperature profiles with different initial helium abundances between
the FG and SG. Our results can provide critical constraints both on
understanding the mixing efficiency in the theoretical models, which is
largely unknown, and on interpretation of the observational carbon abundance
evolution of the bright halo RGB stars.

Stellar populations (1622); Population II stars (1284); Hertzsprung Russell
diagram (725); Globular star clusters (656); Chemical abundances (224);
Stellar evolution (1599); Red giant branch (1368)

††facilities: SMARTS: 1.0 m (STA), WIYN: 0.9 m (HDI, S2KB), Gaia.

##  1 INTRODUCTION

In a generally accepted globular cluster (GC) formation scenario, the second
generation (SG) of the stars formed out of interstellar media polluted by the
first generation (FG) of the stars (e.g., D’ercole et al., 2008). The SG stars
in normal GCs exhibit different elemental abundance patterns than the FG stars
do (Cassisi & Maurizio, 2020; Gratton et al., 2019, and references therein).
For example, the nitrogen enhancement and carbon depletion in the main
sequence (MS) and the lower red giant branch (RGB) stars of the typical GC SGs
could be explained as a natural consequence of the CN cycle occurred in the
previous generation of stars. Photometrically, variations in the carbon and
nitrogen abundances leave distinctive hallmarks through strong absorption band
features of NH, CN, CH molecules, which help to lucidly define multiple
stellar populations (MSPs) of GCs (Lee, 2017, 2019; Lee & Sneden, 2021; Milone
et al., 2017).

Figure 1: (Top) CMDs of M22 in our science field. (Bottom) CMDs for M22
proper-motion member stars.

However, the interpretation of carbon and nitrogen abundances of the bright
RGB stars can be somewhat complicated due to internal mixing episodes
accompanied by their own in situ CN cycle in the hydrogen shell burning
region, which significantly alters their initial carbon and nitrogen
abundances. Many candidates for this non-canonical extra mixing have been
suggested and the thermohaline mixing is considered to be the most promising
one since it appears to reproduce observational constraints, although the
detailed mixing efficiency used in the theoretical models appears to be
largely unknown (e.g., Charbonnel & Zahn, 2007; Lee, 2010; Angelou et al.,
2011).

A few elements, such as 7Li, 12C, 12C/13C, and 14N, have been frequently used
to fine-tune the unknown mixing parameters. It is believed that neither the
[N/Fe] nor 12C/13C are proper probes to investigate the internal mixing, since
the [N/Fe] of the initially nitrogen enhanced stars (i.e., the SG stars) may
not be significantly affected by any extra mixing (Angelou et al., 2011), and
the 12C/13C ratio can rapidly attain near-equilibrium value and quickly
saturated by only moderate amount of mixing (e.g., Sneden, Pilachowski, &
VandenBerg, 1986). The reliable lithium abundance can be obtained with the Li
I resonance doublet at 6707.78 Å in a rather clean spectral region. However,
lithium may not be ideal element to investigate the mixing. Due to its
fragility, lithium can be heavily destructed at the RGB bump (RGBB) luminosity
level and does not provide any useful information for bright GC RGB stars
(e.g., see Figure 1 of Angelou et al., 2015, and references therein). For the
same reason, the degree of lithium depletion is known to be extremely
sensitive to the stellar models (Lattanzio et al., 2015).

The last element standing, carbon, is abundant enough that it cannot be
completely exhausted via the CN cycle in GC RGB stars. Unfortunately, there is
no observable atomic transitions in the optical passband and one should rely
on the diatomic molecular absorption bands, such as CH and CN, to derive
[C/Fe]. A great deal of effort has been directed to determine carbon
abundances by others (e.g., Briley et al., 2004; Smith & Martell, 2003;
Sneden, Pilachowski, & VandenBerg, 1986).

During the past decade, we developed a new photometric system combined with
the robust and self-consistent theoretical fine model grids with various
parameter sets to simultaneously measure the key elements of GC MSPs, [Fe/H],
[C/Fe], and [N/Fe], even in the very crowded field, where traditional
spectroscopic observations cannot be applied (Lee, 2017, 2019, 2021, 2022,
2023; Lee & Sneden, 2021).

In this Letter, we present a photometric [Fe/H], [C/Fe], and [N/Fe] study for
the metal-complex GC M22 (Lee et al., 2009; Lee, 2015, 2016, 2020; Marino et
al., 2009, 2011). With our new [C/Fe] measurements of M22 and those of our
previous studies for M3, M5, and M92, we will discuss the surface carbon
depletion rates of the Milky Way GCs to shed more light on the quantitative
details of the internal mixing processes.

##  2 OBSERVATIONS AND DATA REDUCTION

The journals of observations for M22 are given in Lee (2015, 2020). In 2019,
we obtained additional photometric data for our JWL34 and Strömgren $y$, $b$
for M22, with the total integration times of 6200, 260, and 650 s,
respectively, in three nights from June 27 to July 8 using the KPNO 0.9 m
telescope equipped with a 2 $\times$ 2k CCD chip, providing a field of view
(FOV) of 21′$\times$ 21′.

The raw data handling was described in detail in our previous works (Lee,
2015, 2017). The photometry of M22 and standard stars were analyzed using
DAOPHOTII, DAOGROW, ALLSTAR and ALLFRAME, and COLLECT-CCDAVE-NEWTRIAL packages
(Stetson, 1987, 1994). Finally, we derived the astrometric solutions for
individual stars using the Gaia Early Data Release 3 (EDR3; Gaia
Collaboration, 2020) and the IRAF IMCOORS package.

In order to select M22 member stars, we made use of the proper-motion study
from the Gaia EDR3. We derived the mean proper-motion values using iterative
sigma-clipping calculations, finding that, in units of milliarcsecond per
year, ($\mu_{\rm RA}\times\cos\delta$, $\mu_{\rm decl.}$) = (9.792, $-$5.611)
with standard deviations along the major axis of the ellipse of 0.756 mas yr-1
and along the minor axis of 0.679 mas yr-1. We emphasize that our mean proper
motion for M22 is in good agreement with previous value derived by the Gaia
Collaboration (e.g., Gaia Collaboration, 2018). We considered stars within
3$\sigma$ from the mean values to be M22 proper-motion member stars.

Figure 2: CMDs of member RGB stars with $-$2.5 mag $\leq$ $V-V_{\rm HB}$
$\leq$ 3.0 mag in M22. In the top panels, we show differential reddening
vectors of individual color indices for $\Delta$$E(B-V)$ = 0.1 mag with red
arrows.

##  3 Photometric Indices and Color–Magnitude Diagrams

Throughout this work, we will use our own photometric indices (see also Lee,
2019; Lee & Sneden, 2021; Lee, 2022), defined as

$\displaystyle{{hk}}_{\mathrm{JWL}}$ $\displaystyle=$
$\displaystyle({\mathrm{Ca}}_{\mathrm{JWL}}-b)-(b-y),$ (1)
$\displaystyle{{cn}}_{\mathrm{JWL}}$ $\displaystyle=$
$\displaystyle{JWL}39-{\mathrm{Ca}}_{\mathrm{JWL}},$ (2)
$\displaystyle{{ch}}_{\mathrm{JWL}}$ $\displaystyle=$
$\displaystyle({JWL}43-b)-(b-y),$ (3) $\displaystyle{{nh}}_{\mathrm{JWL}}$
$\displaystyle=$ $\displaystyle({JWL}34-b)-(b-y).$ (4)

The $hk_{\rm JWL}$ index is a good photometric measure of metallicity, while
the $nh_{\rm JWL}$, $cn_{\rm JWL}$, and $ch_{\rm JWL}$ indices are measures of
NH absorption band at $\lambda$3360, CN at $\lambda$3883, and CH at
$\lambda$4250 Å, respectively (e.g., see Lee et al., 2009; Lee, 2015, 2022,
2023; Lee & Sneden, 2021).

In Figure 1, we show our color–magnitude diagrams (CMDs) of M22 using our
color indices. The figure shows that our procedure of selecting cluster’s
member stars by using the Gaia proper motion works excellently and the most of
the foreground and background off-cluster field stars are removed. As we
showed in our previous works (Lee et al., 2009; Lee, 2015, 2020), a bimodal
$hk_{\rm JWL}$ distribution of M22 RGB stars can be clearly seen. Also M22 has
very broad $cn_{\rm JWL}$, $ch_{\rm JWL}$, and $nh_{\rm JWL}$ RGB sequences
due to the variations in the carbon and nitrogen abundances as will be
discussed below (see also Hesser & Harris, 1979; Lee, 2015).

We select the M22 member RGB stars of our interest with $-$2.5 mag $\leq$
$V-V_{\rm HB}$ $\leq$ 3.0 mag in order to derive their [Fe/H]hk, [C/Fe]ch, and
[N/Fe]nh. In Figure 2, we show CMDs of the metal-poor (MP) and the metal-rich
(MR) RGB stars in M22. The classification of the two populations will be
discussed below (see also Lee et al., 2009; Lee, 2015). At a glance, the
differences in the $cn_{\rm JWL}$, and $ch_{\rm JWL}$ CMDs between the MP and
MR populations are noticeable, suggesting that they have different carbon and
nitrogen abundances.

The mean interstellar reddening value toward M22 is rather large, $E(B-V)$ =
0.32 (Harris, 1996, 2010 version), and our M22 CMDs may be vulnerable to
differential reddening effect across our science field. We attempt to correct
the potential differential reddening effect in the following way. First, we
calculated the extinction law for our filter system. We calculated the color
excess of individual color indices using our filter transmission functions,
synthetic spectra of RGB stars for [Fe/H] = $-$1.70 dex with a primordial CNO
and helium abundances as described below, and the interstellar extinction law
by Mathis (1990). Through our calculations, we obtained $E(b-y)$ =
0.820$E(B-V)$, $E$($cn_{\rm JWL}$) = 0.053$E(B-V)$, $E$($ch_{\rm JWL}$) =
$-$0.043$E(B-V)$, $E$($nh_{\rm JWL}$) = 0.539$E(B-V)$, and $E$($hk_{\rm JWL}$)
= $-$0.067$E(B-V)$. In Figure 2, we show differential reddening vectors of
individual color indices for $\Delta$$E(B-V)$ = 0.1 mag. Except for the
$(b-y)$, the reddening vectors do not point to broaden our observed color
indices (i.e., the reddening vectors align along the RGB sequences).
Therefore, it is thought that differential reddening across M22 does not
significantly affect our results (see also Lee, 2015). In order to examine our
claim, we derived the mean $(b-y)$ fiducial sequences of each population and
calculated $\Delta$$E(B-V)$ values of individual RGB stars from their
$E(b-y)$. Finally, we calculated the reddening corrected color indices and we
show them in Figure 2. Except for the $(b-y)$ color, our solutions hardly
reduce the RGB widths, suggesting that the broad RGB widths in $hk_{\rm JWL}$,
$cn_{\rm JWL}$, $ch_{\rm JWL}$, and $nh_{\rm JWL}$ are mainly due to
variations in the metallicity, carbon, and nitrogen abundances.

Figure 3: (a) The [Fe/H]hk distribution of M22 RGB stars. The blue and red
colors denote the MP and MR populations. (b) The histogram of the [Fe/H]hk
along with results returned from our EM estimator. (c) CRDs of each RGB
population. (d–f) Same as (a-c) but for the [Fe/H]hk from the dereddened
$hk_{\rm JWL}$ colors.  Figure 4: (Top) Plots of [C/Fe]ch and [N/Fe]nh of
individual M22 RGB stars against $V$ magnitude. The blue and red colors denote
the MP and MR populations. The error bars in (a) and (b) are for typical
measurement uncertainties in our [C/Fe]ch and [N/Fe]nh derivation,
respectively. (Bottom) Same as the top panels but the results from the
dereddened color indices.

##  4 Metallicity, Carbon, and Nitrogen Abundances

We derive metallicity of individual RGB stars in M22 using our $hk_{\rm JWL}$
measurements. We retrieved the model isochrones for [Fe/H] = $-$2.1, $-$1.9,
$-$1.7, and $-$1.5 dex, $Y$ = 0.247(0.248), 0.275, and 0.300 with
[$\alpha$/Fe] = +0.4 dex, and the age of 12.5 Gyr from a Bag of Stellar Tracks
and Isochrones (Pietrinferni et al., 2021). We adopted different CNO
abundances, [C/Fe] = ($-$0.6, $\Delta$[C/Fe] = 0.2 ,0.6), [N/Fe] = ($-$0.8,
$\Delta$[N/Fe] = 0.4 , 1.6), and [O/Fe] = (0.1, 0.3, 0.5) for each model grid.
Note that our presumed CNO abundances do not affect our photometric
metallicity (Lee, 2022). We constructed 97 model atmospheres and synthetic
spectra for each chemical composition from the lower main sequences to the tip
of RGB sequences using ATLAS12 (Kurucz, 2011) and the latest version of
MOOGSCAT (Sobeck et al., 2011; Sneden, 1974). As we discussed in our previous
works (Lee & Sneden, 2021; Lee, 2023), the latest version of MOOGSCAT (Sobeck
et al., 2011) takes proper care of Rayleigh scattering from neutral hydrogen
(RSNH) atoms with a nonlocal thermodynamic equilibrium treatment of the source
function, which is important to calculate continuum opacities for our short-
wavelength indices such as our JWL34 (i.e., $nh_{\rm JWL}$), due to a
$\lambda^{-4}$ dependency of the RSNH cross section (see Lee & Sneden, 2021,
and references therein).

The photometric metallicity of individual RGB stars can be calculated using
the following relation (also see Appendices of Lee & Sneden, 2021)

$\displaystyle{\rm[Fe/H]}_{hk}$ $\displaystyle\approx$ $\displaystyle
f_{1}({{hk}}_{\mathrm{JWL}},~{}M_{V}).$ (5)

We obtained the mean [Fe/H]hk = $-$1.839 $\pm$ 0.003 dex ($\sigma$ = 0.129)
and we show our results in Figure 3. We emphasize that our photometric
metallicity does not show any gradient against $V$ magnitude and exhibits a
bimodal distribution as already well known (Lee et al., 2009; Lee, 2015, 2016;
Marino et al., 2009, 2011).

In order to perform populational tagging for metallicity, we applied an
expectation-maximization (EM) algorithm for a two-component Gaussian mixture
model on our [Fe/H]hk distribution. Stars with $P$([Fe/H]hk$|x_{i})\geq$ 0.5
from the EM estimator correspond to the MP population, where $x_{i}$ denotes
the individual RGB stars, while those with $P$([Fe/H]hk$|x_{i})$ $<$ 0.5
correspond to the MR population.111Note that the MP of current study is
corresponding to the Ca-w (Lee, 2015) and the G1 (Lee, 2020), while the MR to
the Ca-s and the G2, respectively. We obtained the populational number ratio
of $n$(MP):$n$(MR) = 68.4:31.6 ($\pm$1.4) and our current result is consistent
with our previous result, $n$(Ca-w):$n$(Ca-s) = 70:30 (Lee, 2015). We show the
metallicity distributions of the two populations in Figure 3(b). We obtained
[Fe/H]hk = $-$1.914 $\pm$ 0.002 ($\sigma$ = 0.070) for the MP population and
$-$1.676 $\pm$ 0.002 ($\sigma$ = 0.050) for the MR population. The difference
between the two metallicity groups, $\Delta$[Fe/H]hk = 0.238 $\pm$ 0.003 dex,
is in excellent agreement with those in spectroscopic analyses (e.g., Lee,
2016). We also examined the cumulative radial distributions (CRDs) and we show
them in Figure 3(c). We performed Kolmogorov–Smirnov tests and obtained that
the MP and MR populations are most likely drawn from the different parent CRDs
(see also Lee et al., 2009; Lee, 2015).

With our differential reddening estimates as mentioned above, we calculate the
metallicity from the dereddened $hk_{\rm JWL}$ colors, finding [Fe/H]hk =
$-$1.840 $\pm$ 0.003 dex ($\sigma$ = 0.129), $-$1.915 $\pm$ 0.002 ($\sigma$ =
0.069), $-$1.676 $\pm$ 0.002 ($\sigma$ = 0.058) for the mean value, the MP,
and MR populations with the identical populational number ratio
$n$(MP):$n$(MR) = 68.4:31.6 ($\pm$1.4). As we pointed out above, we emphasize
again that the differential reddening does not appear to affect our results.

Using the following relations (Lee, 2021; Lee & Sneden, 2021; Lee, 2023),

$\displaystyle{\rm[C/Fe]}_{ch}$ $\displaystyle\approx$ $\displaystyle
f_{2}({{ch}}_{\mathrm{JWL}},~{}{\mathrm{[Fe/H]}_{hk}},~{}M_{V}),$ (6)
$\displaystyle{\rm[N/Fe]}_{nh}$ $\displaystyle\approx$ $\displaystyle
f_{3}({{nh}}_{\mathrm{JWL}},~{}{\mathrm{[Fe/H]}_{hk}},~{}M_{V}),$ (7)

we derive the photometric [C/Fe]ch and [N/Fe]nh of each population in M22. In
Figure 4, we show [C/Fe]ch and [N/Fe]nh against the $V$ magnitude with and
without the differential reddening correction, suggesting that the
differential reddening correction does not significantly affect our
photometric [C/Fe]ch and [N/Fe]nh. Our results clearly show that the carbon
abundance decreases and nitrogen abundance increases in RGB stars brighter
than the RGBB ($V$ $\approx$ 14.0 mag for M22; Lee 2015) as can be seen in our
previous studies of other Galactic GCs, M5, M3, and M92 (Lee & Sneden, 2021;
Lee, 2021, 2023).

In Figures 4(i–j), the distributions of the MP and MR on the [C/Fe]ch versus
[N/Fe]nh plane are different in the sense that at a given carbon abundance the
nitrogen abundance of the MP RGB stars are lower, confirming previous results
of Marino et al. (2011) and Lee (2015). As Marino et al. (2011) noted, without
any discernible oxygen abundance differences, the difference in the total CNO
abundances between the two populations in M22 is likely responsible for this
separation.

In Figure 4, we emphasize that the degrees of the carbon depletion between the
M22 MP and MR RGB stars brighter than the RGBB appears to be quite different.
On the other hand, the degree of the nitrogen enhancement is not as clear as
that of the carbon depletion, at least due to the dependence of the nitrogen
enhancement on their initial nitrogen abundances (e.g., Angelou et al., 2011).

We performed subpopulational tagging for the MP and MR populations. Due to the
carbon depletion and the nitrogen enhancement in RGB stars brighter than the
RGBB, neither [C/Fe]ch nor [N/Fe]nh may be the proper probes for the
populational tagging. Instead, we perform populational tagging using our
$\parallel$$ch_{\rm JWL}$ and $\parallel$$nh_{\rm JWL}$, which are defined as

$\displaystyle\parallel ch_{\rm JWL}$ $\displaystyle\equiv$
$\displaystyle\frac{ch_{\rm JWL}-ch_{\rm JWL,red}}{ch_{\rm JWL,red}-ch_{\rm
JWL,blue}},$ (8) $\displaystyle\parallel nh_{\rm JWL}$ $\displaystyle\equiv$
$\displaystyle\frac{nh_{\rm JWL}-nh_{\rm JWL,red}}{nh_{\rm JWL,red}-nh_{\rm
JWL,blue}},$ (9)

where the subscripts denote the fiducials of the red and the blue sequences of
individual color indices (see also Lee & Sneden, 2021). With these procedures,
the curvatures of the RGB sequence on the $ch_{\rm JWL}$ and $nh_{\rm JWL}$
CMDs can be effectively removed. Then we can perform populational tagging for
the bright RGB stars on these normalized color indices (e.g., Lee, 2020). We
made histograms for the ($\parallel$$ch_{\rm JWL}$ $-$ $\parallel$$nh_{\rm
JWL}$), which is roughly equivalent to [C/N], showing at least two distinctive
peaks in each histogram. We applied EM algorithms for each population and
obtained subpopulational number ratios of $n$(FG):$n$(SG) = 49.0:51.0
($\pm$2.4) for the MP and 50.0:50.0 ($\pm$3.5) for the MR.

With our subpopulations, we derived $V$ magnitude dependence of carbon
abundance depletion, $d\mathrm{[C/Fe]}/M_{V}$. In order to derive the slopes
of such relation, we used a ordinary least-squares fit and a robust fit that
minimizes absolute deviation (MEDFIT: Press et al., 1986). We show our results
in Figure 5 and Table 1. As shown, the M22 MP population has steeper gradients
than the MR does, mainly due to difference in the mean metallicity between the
two populations (Charbonnel & Zahn, 2007). Also importantly, the SG
subpopulations appear to have slightly steeper gradients than the FG, although
the differences are not statistically significant.

Figure 5: Plots of [C/Fe]ch vs. $M_{V}$ of individual populations of M22, M5, M3, and M92 for RGB stars brighter than their RGBB $V$ magnitudes. The blue and red solid lines denote the least-squares fitting and the robust fitting by minimizing absolute deviation, respectively. The numbers in parentheses are the mean metallicity of each GC. Table 1: Magnitude Dependence of Carbon Abundance. | Pop. | [Fe/H]aaHarris (1996, 2010 version) | [Fe/H]hk |  $d$[C/Fe]/$dM_{V}$bbA least-squares fitting. |  $d$[C/Fe]/$dM_{V}$ccA fitting by minimizing absolute deviation. |  $d$[C/Fe]/$dM_{V}$ddSmith & Martell (2003).
---|---|---|---|---|---|---
M22-MR | FG |  |  $-$1.68 | 0.104 $\pm$ 0.033 | 0.150 $\pm$ 0.088 |
| SG |  |  | 0.142 $\pm$ 0.029 | 0.156 $\pm$ 0.074 |
M22-MP | FG |  |  $-$1.91 | 0.203 $\pm$ 0.024 | 0.235 $\pm$ 0.090 |
| SG |  |  | 0.227 $\pm$ 0.020 | 0.198 $\pm$ 0.092 |
M5 | FG |  $-$1.29 |  $-$1.30 | 0.096 $\pm$ 0.019 | 0.077 $\pm$ 0.056 |
| SG |  |  | 0.106 $\pm$ 0.017 | 0.141 $\pm$ 0.089 |
M3 | FG |  $-$1.50 |  $-$1.50 | 0.159 $\pm$ 0.012 | 0.105 $\pm$ 0.049 | 0.236 $\pm$ 0.033
| SG |  |  | 0.206 $\pm$ 0.016 | 0.187 $\pm$ 0.055 |
M92 | FG |  $-$2.31 |  $-$2.31 | 0.391 $\pm$ 0.022 | 0.315 $\pm$ 0.066 | 0.227 $\pm$ 0.045
| SG |  |  | 0.374 $\pm$ 0.032 | 0.396 $\pm$ 0.119 |
Figure 6: (a) A plot of $d\mathrm{[C/Fe]}/M_{V}$ vs. [Fe/H] returned from the
least-squares fitting (blue solid lines in Figure 5). The blue and red colors
denote the FG and SG, respectively. The green and gold solid lines show a
linear regression with confidence intervals of 95%. The yellow shade box
indicates the mean value by Smith & Martell (2003). (b) Same as (a) but for
results returned from the robust fitting by minimizing absolute deviation (red
solid lines in Figure 5).

##  5 Surface Carbon Depletion Rates of Milky Way Globular Clusters

It is well known that the mixing efficiency depends on metallicity since the
hydrogen burning shell will be hotter and wider with decreasing metallicity
(e.g., Church et al., 2014). We investigate the degree of the carbon depletion
against $V$ magnitude, $d\mathrm{[C/Fe]}/M_{V}$, with expand samples of our
previous studies, M5 ([Fe/H] = $-$1.29; Lee, 2021), M3 ($-$1.50; Lee & Sneden,
2021), and M92 ($-$2.31; Lee, 2023). As with M22, we derive the slopes of the
relations using two different methods. We show the results in Figure 5 and
Table 1. In the table, the $d\mathrm{[C/Fe]}/M_{V}$ of M3 and M92 by Smith &
Martell (2003) are also listed. Smith & Martell (2003) did not distinguish the
FG and SG in their analyses and their slopes for M3 and M92 are slightly
different from our values.

In Figure 6, we show plots of $d\mathrm{[C/Fe]}/M_{V}$ versus [Fe/H] for the
case of a least-squares fit (Case 1; column (5) of Table 1) and a robust fit,
which minimizes absolute deviation (Case 2; column (6) of Table 1),
respectively. We obtained following relations

$\displaystyle\frac{d\mathrm{[C/Fe]}}{dM_{V}}$ $\displaystyle=$
$\displaystyle-0.261(\pm 0.043)\mathrm{[Fe/H]}-0.253(\pm 0.077),$ (10)
$\displaystyle\frac{d\mathrm{[C/Fe]}}{dM_{V}}$ $\displaystyle=$
$\displaystyle-0.241(\pm 0.037)\mathrm{[Fe/H]}-0.224(\pm 0.066),$ (11)

and, not surprisingly, our results clearly show that the
$d\mathrm{[C/Fe]}/M_{V}$ depends on metallicity. Note that Smith & Martell
(2003) obtained the mean value of $d\mathrm{[C/Fe]}/M_{V}$ $\approx$ 0.22
$\pm$ 0.03 from three GCs (M92, NGC 6397, and M3) and the MP halo giants. On
the other hand, Martell, Smith, & Briley (2008) argued that the carbon
depletion rate222We note that the carbon depletion rate by Martell, Smith, &
Briley (2008), $\Delta$[C/Fe]/$\Delta$t, is different from our value,
$d\mathrm{[C/Fe]}/M_{V}$. becomes doubled from [Fe/H] = $-$1.3 to $-$2.3.

As we mentioned above, the slopes of the $d\mathrm{[C/Fe]}/M_{V}$ for the FG
and SG appear to be slightly different for M5, M3, and M92, in the sense that
the SG tends to have a steeper slope. We argue that different initial helium
abundances between the FG and SG are responsible for the difference in the
slopes. Due to the high helium abundances of the SG stars, which can be
inferred from their RGBB $V$ magnitudes (e.g., Lee, 2020, 2023; Lee & Sneden,
2021), the SG RGB stars have high temperature at the hydrogen burning shell,
which causes slightly faster carbon destruction through the CN cycle (e.g.,
Church et al., 2014).

##  6 SUMMARY

We investigated the photometric metallicity, carbon, and nitrogen abundances
of the metal-complex globular cluster M22. Our results confirmed previous
results that M22 contains at least two MSPs with heterogeneous metallicities
(Lee et al., 2009; Lee, 2015, 2016; Marino et al., 2009, 2011).

We obtained [Fe/H]hk = $-$1.839 $\pm$ 0.003 dex ($\sigma$ = 0.129), $-$1.914
$\pm$ 0.002 ($\sigma$ = 0.070), and $-$1.676 $\pm$ 0.002 ($\sigma$ = 0.050)
for the all RGB, the MP, and MR RGB populations, respectively. Applying
differential reddening correction does not affect our photometric metallicity
measurements, no larger than 0.001 dex in the mean values.

Our [C/Fe]ch, [N/Fe]nh measurements of individual populations in M22 showed
that each population showed evidences of the CN process accompanied by deep
mixing episodes with different degrees.

With our carbon abundance measurements for M22 and other GCs in our previous
studies, M5 (Lee, 2021), M3 (Lee & Sneden, 2021), M92 (Lee, 2023), we
investigated the surface carbon depletion rates, $d\mathrm{[C/Fe]}/M_{V}$, for
Milky Way GCs against metallicity, finding $d\mathrm{[C/Fe]}/M_{V}$ $\propto$
$-$0.25[Fe/H]. We also argued that the carbon depletion rates of the SG are
larger than those of the FG, most likely due to different initial helium
abundances between the FG and SG, which cause different internal temperature
profiles.

Our results presented here can provide critical constraints both on
understanding the mixing efficiency in the theoretical models, which is
largely unknown, and on interpretation of the carbon abundances evolution seen
from the bright halo RGB stars.

J.-W.L. acknowledges financial support from the Basic Science Research Program
(grant No. 2019R1A2C2086290) through the National Research Foundation of Korea
(NRF) and from the faculty research fund of Sejong University in 2022. He also
thanks the anonymous referee for encouraging comments.

## References

  * Angelou et al. (2011) Angelou, G. C., Church, R. P., STancliffe, R. J., Lattanzio, J. C., & Smith, G. H. 2011, ApJ, 728, 79
  * Angelou et al. (2015) Angelou, G. C., D’Orazi, V., Constantino, T. N., Church, R. P., et al. 2015, MNRAS, 450, 2423
  * Briley et al. (2004) Briley, M. M., Harbeck, D., Smith, G. H., & Grebel, E. K. AJ, 127, 1588
  * Cassisi & Maurizio (2020) Cassisi, S. & Maurizio, S. 2020, A&A Rev., 28, 5
  * Charbonnel & Zahn (2007) Charbonnel, C., & Zahn, J.-P. 2007, A&A, 467, L15
  * Church et al. (2014) Church, R. P., Lattanzio, J., Angelou, G., Tout, C. A., & STancliffe, R. J. 2014, MNRAS, 443, 977
  * D’ercole et al. (2008) D’Ercole, A., Vesperini, E., D’Antona, F., McMillan, S. L. W., & Recchi, S. 2008, MNRAS, 391, 825
  * Gaia Collaboration (2020) Gaia Collaboration, Brown, A. G. A., Vallenari, A., et al. 2021, A&A, 649, 1
  * Gaia Collaboration (2018) Gaia Collaboration, Helmi, A., et al. 2018, A&A, 616, 12
  * Gratton et al. (2019) Gratton, R., Bragaglia, A., Carretta, E., D’Orazi, V. et al. 2019, A&A Rev., 27, 8
  * Harris (1996) Harris, W. E. 1996, AJ, 112, 1487
  * Hesser & Harris (1979) Hesser, J .E., & Harris, W. E. 1979, ApJ, 234, 513
  * Kurucz (2011) Kurucz, R. L. 2011, Can. J. Phys., 89, 417
  * Lattanzio et al. (2015) Lattanzio, J. C., Siess, L., Church, R. P., Angelou, G., et al. 2015, MNRAS, 446, 2673
  * Lee (2010) Lee, J.-W. 2010, MNRAS, 405, L36
  * Lee (2015) Lee, J.-W. 2015, ApJS, 219, 7
  * Lee (2016) Lee, J.-W. 2016, ApJS, 226, 16
  * Lee (2017) Lee, J.-W. 2017, ApJ, 844, 77
  * Lee (2019) Lee, J.-W. 2019, ApJ, 872, 41
  * Lee (2020) Lee, J.-W. 2020, ApJ, 888, L6
  * Lee (2021) Lee, J.-W. 2021, ApJ, 918, L24
  * Lee (2022) Lee, J.-W. 2022, ApJS, 263, 20
  * Lee (2023) Lee, J.-W. 2023, ApJ, 948, L16
  * Lee et al. (2009) Lee, J.-W., Kang, Y.-W., Lee, J., & Lee, Y.-W. 2009, Nature, 462, 480
  * Lee & Sneden (2021) Lee, J.-W., & Sneden, C. 2021, ApJ, 909, 167
  * Marino et al. (2011) Marino, A. F., Sneden, C., Kraft, R. P., et al. 2011, A&A, 532, A8
  * Marino et al. (2009) Marino, A. F., Milone, A. P., Piotto, G., et al. 2009, A&A, 505, 1099
  * Martell, Smith, & Briley (2008) Martell, S. L., Smith, G. H., Briely, M. 2009, AJ, 136, 2522
  * Mathis (1990) Mathis, J. S. 1999 ARA&A, 28, 37
  * Milone et al. (2017) Milone A. P., Piotto, G., Renzini, A. et al. 2017, MNRAS, 464, 3636
  * Pietrinferni et al. (2021) Pietrinferni, A., Hidalgo, S., Cassisi, S., Salaris, M., et al. 2021, ApJ, 908, 102
  * Press et al. (1986) Press, W. H., Teukolsky, S. A., Vetterling, W. T., & Flannery, B. P. 1992, Numerical Recipes in Fortran 77 : The Art of Scientific Computing (Cambridge: Cambridge University Press)
  * Smith & Martell (2003) Smith, G. H., & Martell, S. L. 2003, PASP, 115, 1211
  * Sneden (1974) Sneden, C. 1974, ApJ, 189, 493
  * Sneden, Pilachowski, & VandenBerg (1986) Sneden, C., Pilachowski, C. A., & VandenBerg, D. A., 1986, ApJ, 311, 816
  * Sobeck et al. (2011) Sobeck, J. E., Kraft, R. P., Sneden, C., et al. 2011, AJ, 141, 175
  * Stetson (1987) Stetson P. B. 1987, PASP, 99, 191
  * Stetson (1994) Stetson P. B. 1994, PASP, 106, 250


††thanks: Deceased††thanks: Deceased††thanks: Deceased††thanks: Deceased

The BABAR Collaboration

# Search for Lepton-Flavor-Violating Decays $D^{0}\rightarrow
X^{0}e^{\pm}\mu^{\mp}$

J. P. Lees  V. Poireau  V. Tisserand  Laboratoire d’Annecy-le-Vieux de
Physique des Particules (LAPP), Université de Savoie, CNRS/IN2P3, F-74941
Annecy-Le-Vieux, France  E. Grauges  Universitat de Barcelona, Facultat de
Fisica, Departament ECM, E-08028 Barcelona, Spain  A. Palano  INFN Sezione di
Bari and Dipartimento di Fisica, Università di Bari, I-70126 Bari, Italy  G.
Eigen  University of Bergen, Institute of Physics, N-5007 Bergen, Norway  D.
N. Brown  Yu. G. Kolomensky  Lawrence Berkeley National Laboratory and
University of California, Berkeley, California 94720, USA  M. Fritsch  H. Koch
T. Schroeder  Ruhr Universität Bochum, Institut für Experimentalphysik 1,
D-44780 Bochum, Germany  R. Cheaibb C. Heartyab T. S. Mattisonb J. A. McKennab
R. Y. Sob Institute of Particle Physics${}^{\,a}$; University of British
Columbiab, Vancouver, British Columbia, Canada V6T 1Z1  V. E. Blinovabc A. R.
Buzykaeva V. P. Druzhininab V. B. Golubevab E. A. Kozyrevab E. A. Kravchenkoab
A. P. Onuchinabc S. I. Serednyakovab Yu. I. Skovpenab E. P. Solodovab K. Yu.
Todyshevab Budker Institute of Nuclear Physics SB RAS, Novosibirsk 630090a,
Novosibirsk State University, Novosibirsk 630090b, Novosibirsk State Technical
University, Novosibirsk 630092c, Russia  A. J. Lankford  University of
California at Irvine, Irvine, California 92697, USA  B. Dey  J. W. Gary  O.
Long  University of California at Riverside, Riverside, California 92521, USA
A. M. Eisner  W. S. Lockman  W. Panduro Vazquez  University of California at
Santa Cruz, Institute for Particle Physics, Santa Cruz, California 95064, USA
D. S. Chao  C. H. Cheng  B. Echenard  K. T. Flood  D. G. Hitlin  J. Kim  Y. Li
D. X. Lin  T. S. Miyashita  P. Ongmongkolkul  J. Oyang  F. C. Porter  M.
Röhrken  California Institute of Technology, Pasadena, California 91125, USA
Z. Huard  B. T. Meadows  B. G. Pushpawela  M. D. Sokoloff  L. Sun  Now at:
Wuhan University, Wuhan 430072, China  University of Cincinnati, Cincinnati,
Ohio 45221, USA  J. G. Smith  S. R. Wagner  University of Colorado, Boulder,
Colorado 80309, USA  D. Bernard  M. Verderi  Laboratoire Leprince-Ringuet,
Ecole Polytechnique, CNRS/IN2P3, F-91128 Palaiseau, France  D. Bettonia C.
Bozzia R. Calabreseab G. Cibinettoab E. Fioravantiab I. Garziaab E. Luppiab V.
Santoroa INFN Sezione di Ferraraa; Dipartimento di Fisica e Scienze della
Terra, Università di Ferrarab, I-44122 Ferrara, Italy  A. Calcaterra  R. de
Sangro  G. Finocchiaro  S. Martellotti  P. Patteri  I. M. Peruzzi  M. Piccolo
M. Rotondo  A. Zallo  INFN Laboratori Nazionali di Frascati, I-00044 Frascati,
Italy  S. Passaggio  C. Patrignani  Now at: Università di Bologna and INFN
Sezione di Bologna, I-47921 Rimini, Italy  INFN Sezione di Genova, I-16146
Genova, Italy  B. J. Shuve  Harvey Mudd College, Claremont, California 91711,
USA  H. M. Lacker  Humboldt-Universität zu Berlin, Institut für Physik,
D-12489 Berlin, Germany  B. Bhuyan  Indian Institute of Technology Guwahati,
Guwahati, Assam, 781 039, India  U. Mallik  University of Iowa, Iowa City,
Iowa 52242, USA  C. Chen  J. Cochran  S. Prell  Iowa State University, Ames,
Iowa 50011, USA  A. V. Gritsan  Johns Hopkins University, Baltimore, Maryland
21218, USA  N. Arnaud  M. Davier  F. Le Diberder  A. M. Lutz  G. Wormser
Université Paris-Saclay, CNRS/IN2P3, IJCLab, F-91405 Orsay, France  D. J.
Lange  D. M. Wright  Lawrence Livermore National Laboratory, Livermore,
California 94550, USA  J. P. Coleman  E. Gabathuler  D. E. Hutchcroft  D. J.
Payne  C. Touramanis  University of Liverpool, Liverpool L69 7ZE, United
Kingdom  A. J. Bevan  F. Di Lodovico  Now at: King’s College, London, WC2R
2LS, UK  R. Sacco  Queen Mary, University of London, London, E1 4NS, United
Kingdom  G. Cowan  University of London, Royal Holloway and Bedford New
College, Egham, Surrey TW20 0EX, United Kingdom  Sw. Banerjee  D. N. Brown  C.
L. Davis  University of Louisville, Louisville, Kentucky 40292, USA  A. G.
Denig  W. Gradl  K. Griessinger  A. Hafner  K. R. Schubert  Johannes
Gutenberg-Universität Mainz, Institut für Kernphysik, D-55099 Mainz, Germany
R. J. Barlow  Now at: University of Huddersfield, Huddersfield HD1 3DH, UK  G.
D. Lafferty  University of Manchester, Manchester M13 9PL, United Kingdom  R.
Cenci  A. Jawahery  D. A. Roberts  University of Maryland, College Park,
Maryland 20742, USA  R. Cowan  Massachusetts Institute of Technology,
Laboratory for Nuclear Science, Cambridge, Massachusetts 02139, USA  S. H.
Robertsonab R. M. Seddonb Institute of Particle Physics${}^{\,a}$; McGill
Universityb, Montréal, Québec, Canada H3A 2T8  N. Neria F. Palomboab INFN
Sezione di Milanoa; Dipartimento di Fisica, Università di Milanob, I-20133
Milano, Italy  L. Cremaldi  R. Godang  Now at: University of South Alabama,
Mobile, Alabama 36688, USA  D. J. Summers  University of Mississippi,
University, Mississippi 38677, USA  P. Taras  Université de Montréal, Physique
des Particules, Montréal, Québec, Canada H3C 3J7  G. De Nardo  C. Sciacca
INFN Sezione di Napoli and Dipartimento di Scienze Fisiche, Università di
Napoli Federico II, I-80126 Napoli, Italy  G. Raven  NIKHEF, National
Institute for Nuclear Physics and High Energy Physics, NL-1009 DB Amsterdam,
The Netherlands  C. P. Jessop  J. M. LoSecco  University of Notre Dame, Notre
Dame, Indiana 46556, USA  K. Honscheid  R. Kass  Ohio State University,
Columbus, Ohio 43210, USA  A. Gaza M. Margoniab M. Posoccoa G. Simiab F.
Simonettoab R. Stroiliab INFN Sezione di Padovaa; Dipartimento di Fisica,
Università di Padovab, I-35131 Padova, Italy  S. Akar  E. Ben-Haim  M. Bomben
G. R. Bonneaud  G. Calderini  J. Chauveau  G. Marchiori  J. Ocariz
Laboratoire de Physique Nucléaire et de Hautes Energies, Sorbonne Université,
Paris Diderot Sorbonne Paris Cité, CNRS/IN2P3, F-75252 Paris, France  M.
Biasiniab E. Manonia A. Rossia INFN Sezione di Perugiaa; Dipartimento di
Fisica, Università di Perugiab, I-06123 Perugia, Italy  G. Batignaniab S.
Bettariniab M. Carpinelliab Also at: Università di Sassari, I-07100 Sassari,
Italy  G. Casarosaab M. Chrzaszcza F. Fortiab M. A. Giorgiab A. Lusianiac B.
Oberhofab E. Paoloniab M. Ramaa G. Rizzoab J. J. Walsha L. Zaniab INFN Sezione
di Pisaa; Dipartimento di Fisica, Università di Pisab; Scuola Normale
Superiore di Pisac, I-56127 Pisa, Italy  A. J. S. Smith  Princeton University,
Princeton, New Jersey 08544, USA  F. Anullia R. Facciniab F. Ferrarottoa F.
Ferronia Also at: Gran Sasso Science Institute, I-67100 L’Aquila, Italy  A.
Pilloniab G. Pireddaa INFN Sezione di Romaa; Dipartimento di Fisica,
Università di Roma La Sapienzab, I-00185 Roma, Italy  C. Bünger  S. Dittrich
O. Grünberg  M. Heß  T. Leddig  C. Voß  R. Waldi  Universität Rostock, D-18051
Rostock, Germany  T. Adye  F. F. Wilson  Rutherford Appleton Laboratory,
Chilton, Didcot, Oxon, OX11 0QX, United Kingdom  S. Emery  G. Vasseur  IRFU,
CEA, Université Paris-Saclay, F-91191 Gif-sur-Yvette, France  D. Aston  C.
Cartaro  M. R. Convery  J. Dorfan  W. Dunwoodie  M. Ebert  R. C. Field  B. G.
Fulsom  M. T. Graham  C. Hast  W. R. Innes  P. Kim  D. W. G. S. Leith  S.
Luitz  D. B. MacFarlane  D. R. Muller  H. Neal  B. N. Ratcliff  A. Roodman  M.
K. Sullivan  J. Va’vra  W. J. Wisniewski  SLAC National Accelerator
Laboratory, Stanford, California 94309 USA  M. V. Purohit  J. R. Wilson
University of South Carolina, Columbia, South Carolina 29208, USA  A. Randle-
Conde  S. J. Sekula  Southern Methodist University, Dallas, Texas 75275, USA
H. Ahmed  St. Francis Xavier University, Antigonish, Nova Scotia, Canada B2G
2W5  M. Bellis  P. R. Burchat  E. M. T. Puccio  Stanford University, Stanford,
California 94305, USA  M. S. Alam  J. A. Ernst  State University of New York,
Albany, New York 12222, USA  R. Gorodeisky  N. Guttman  D. R. Peimer  A.
Soffer  Tel Aviv University, School of Physics and Astronomy, Tel Aviv, 69978,
Israel  S. M. Spanier  University of Tennessee, Knoxville, Tennessee 37996,
USA  J. L. Ritchie  R. F. Schwitters  University of Texas at Austin, Austin,
Texas 78712, USA  J. M. Izen  X. C. Lou  University of Texas at Dallas,
Richardson, Texas 75083, USA  F. Bianchiab F. De Moriab A. Filippia D. Gambaab
INFN Sezione di Torinoa; Dipartimento di Fisica, Università di Torinob,
I-10125 Torino, Italy  L. Lanceri  L. Vitale  INFN Sezione di Trieste and
Dipartimento di Fisica, Università di Trieste, I-34127 Trieste, Italy  F.
Martinez-Vidal  A. Oyanguren  IFIC, Universitat de Valencia-CSIC, E-46071
Valencia, Spain  J. Albertb A. Beaulieub F. U. Bernlochnerb G. J. Kingb R.
Kowalewskib T. Lueckb I. M. Nugentb J. M. Roneyb R. J. Sobieab N. Tasneemb
Institute of Particle Physics${}^{\,a}$; University of Victoriab, Victoria,
British Columbia, Canada V8W 3P6  T. J. Gershon  P. F. Harrison  T. E. Latham
Department of Physics, University of Warwick, Coventry CV4 7AL, United Kingdom
R. Prepost  S. L. Wu  University of Wisconsin, Madison, Wisconsin 53706, USA

###### Abstract

We present a search for seven lepton-flavor-violating neutral charm meson
decays of the type $D^{0}\rightarrow X^{0}e^{\pm}\mu^{\mp}$, where $X^{0}$
represents a $\pi^{0}$, $K^{0}_{\rm S}$, ${{\kern
1.99997pt\overline{\kern-1.99997ptK}{}}}{}^{*0}$, $\rho^{0}$, $\phi$,
$\omega$, or $\eta$ meson. The analysis is based on $468$ $\mbox{\,fb}^{-1}$
of $e^{+}e^{-}$ annihilation data collected at or close to the
${\it\Upsilon}(4S)$ resonance with the BABAR detector at the SLAC National
Accelerator Laboratory. No significant signals are observed, and we establish
90% confidence level upper limits on the branching fractions in the range
$(5.0-22.5)\times 10^{-7}$. The limits are between 1 and 2 orders of magnitude
more stringent than previous measurements.

###### pacs:

13.25.Ft, 11.30.Fs

††preprint: BABAR-PUB-20/001††preprint: SLAC-PUB-17524

##  I Introduction

Lepton-flavor-conserving charm decays such as $D\rightarrow X{e^{+}e^{-}}$ or
$D\rightarrow X{\mu^{+}\mu^{-}}$, where $X$ is a meson, can occur in the
standard model (SM) through short-distance Paul  _et al._ (2011); Schwartz
(1993) and long-distance Schwartz (1993) processes, which have branching
fractions of order ${\mathcal{O}}(10^{-9})$ and ${\mathcal{O}}(10^{-6})$,
respectively. In contrast, the lepton-flavor-violating (LFV) neutral charm
decays ${{D}^{0}}\rightarrow X^{0}e^{\pm}\mu^{\mp}$, where $X^{0}$ is a
neutral meson, are effectively forbidden in the SM because they can occur only
through lepton-flavor mixing Guadagnoli and Lane (2015) and are therefore
suppressed to the order ${\mathcal{O}}(10^{-50})$. As such, the decays
${{D}^{0}}\rightarrow X^{0}e^{\pm}\mu^{\mp}$ should not be visible with
current data samples. However, new-physics models, such as those involving
Majorana neutrinos, leptoquarks, and two-Higgs doublets, allow for lepton
number and lepton flavor to be violated de Boer and Hiller (2016); Fajfer and
Košnik (2015); Atre  _et al._ (2009); Yuan  _et al._ (2013); Hai-Rong  _et
al._ (2015). Some models make predictions for, or use constraints from, three-
body decays of the form ${D}\rightarrow Xl^{\prime}l$ or ${B}\rightarrow
Xl^{\prime}l$, where $\ell$ and $\ell^{\prime}$ represent an electron or muon
Paul  _et al._ (2011, 2014); Burdman  _et al._ (2002); Fajfer and Prelovšek
(2006); Fajfer  _et al._ (2007); Atre  _et al._ (2009); Yuan  _et al._ (2013).
Most recent theoretical work has targeted the charged charm decays
${{D}^{+}}\rightarrow X^{+}\ell^{\prime+}{\ell^{-}}$. For example, Ref. de
Boer and Hiller (2016) estimates that
${\mathcal{B}}({{D}^{+}}\rightarrow{{\pi}^{+}}\mu^{\pm}e^{\mp})$ can be as
large as $2\times 10^{-6}$ for certain leptoquark couplings. Some models that
consider LFV and lepton-number-violating (LNV) four-body charm decays, with
two leptons and two hadrons in the final state, predict branching fractions up
to ${\mathcal{O}}(10^{-5})$, approaching those accessible with current data
Atre  _et al._ (2009); Yuan  _et al._ (2013); Hai-Rong  _et al._ (2015).

The branching fractions ${\mathcal{B}}({{D}^{0}}\rightarrow
h^{\prime-}h^{+}{\mu^{+}\mu^{-}})$, where $h^{\prime}$ and $h$ represent a $K$
or $\pi$ meson, and
${\mathcal{B}}({{D}^{0}}\rightarrow{{K}^{-}}{{\pi}^{+}}{e^{+}e^{-}})$ have
recently been measured to be ${\mathcal{O}}(10^{-7})$ to
${\mathcal{O}}(10^{-6})$ Aaij  _et al._ (2016, 2017); Lees  _et al._ (2019),
compatible with SM predictions Cappiello  _et al._ (2013); de Boer and Hiller
(2018). The branching fractions for the decays ${{D}^{0}}\rightarrow
X^{0}{e^{+}e^{-}}$ and ${{D}^{0}}\rightarrow X^{0}{\mu^{+}\mu^{-}}$ have not
yet been measured. However 90% confidence level (C.L.) upper limits on the
branching fractions do exist and are in the range $(0.3-10)\times 10^{-5}$ for
${{D}^{0}}\rightarrow X^{0}{e^{+}e^{-}}$ and $(3.2-53)\times 10^{-5}$ for
${{D}^{0}}\rightarrow X^{0}{\mu^{+}\mu^{-}}$ Kodama  _et al._ (1995);
Freyberger  _et al._ (1996); Aitala  _et al._ (2001); Ablikim  _et al._
(2018). It is likely that one or more of these decays are a major contributor
to the branching fractions of the decays ${{D}^{0}}\rightarrow
h^{\prime-}h^{+}{e^{+}e^{-}}$ or ${{D}^{0}}\rightarrow
h^{\prime-}h^{+}{\mu^{+}\mu^{-}}$, as long-distance processes are predicted to
be dominant Schwartz (1993), and published distributions of the invariant
masses $m(h^{\prime-}h^{+})$ for ${{D}^{0}}\rightarrow
h^{\prime-}h^{+}{\mu^{+}\mu^{-}}$ and
${{D}^{0}}\rightarrow{{K}^{-}}{{\pi}^{+}}{e^{+}e^{-}}$ indicate large yields
near some of the $X^{0}$ invariant masses Aaij  _et al._ (2016, 2017); Lees
_et al._ (2019).

The most stringent existing upper limits on the branching fractions for the
LFV four-body decays of the type ${{D}^{0}}\rightarrow
h^{\prime-}h^{+}e^{\pm}\mu^{\mp}$ are in the range $(11.0-19.0)\times 10^{-7}$
at the 90% confidence level Lees  _et al._ (2020). For the LFV decays
${{D}^{0}}\rightarrow X^{0}e^{\pm}\mu^{\mp}$, where $X^{0}$ is an intermediate
resonance meson decaying to $h^{\prime-}h^{+}$,
${{\pi}^{+}}{{\pi}^{-}}{{\pi}^{0}}$ or $\gamma\gamma$, the 90% C.L. limits are
in the range $(3.4-118)\times 10^{-5}$ Freyberger  _et al._ (1996); Aitala
_et al._ (2001); Tanabashi  _et al._ (2018). For the ${{D}^{0}}\rightarrow
X^{0}e^{\pm}\mu^{\mp}$ decays with the same final state as the
${{D}^{0}}\rightarrow h^{\prime-}h^{+}e^{\pm}\mu^{\mp}$ decays
(${{D}^{0}}\rightarrow{{K}^{0}_{\mathrm{\scriptscriptstyle
S}}}(\rightarrow{{\pi}^{+}}{{\pi}^{-}})e^{\pm}\mu^{\mp}$,
${{D}^{0}}\rightarrow{{\rho}^{0}}(\rightarrow{{\pi}^{+}}{{\pi}^{-}})e^{\pm}\mu^{\mp}$,
${{D}^{0}}\rightarrow{{\kern
1.99997pt\overline{\kern-1.99997ptK}{}}{}^{*0}}(\rightarrow{{K}^{-}}{{\pi}^{+}})e^{\pm}\mu^{\mp}$,
and
${{D}^{0}}\rightarrow\phi(\rightarrow{{K}^{+}}{{K}^{-}})e^{\pm}\mu^{\mp}$),
the current ${{D}^{0}}\rightarrow X^{0}e^{\pm}\mu^{\mp}$ branching fraction
upper limits, which are in the range $(3.4-8.3)\times 10^{-5}$ Freyberger  _et
al._ (1996); Aitala  _et al._ (2001); Tanabashi  _et al._ (2018), are
approximately 20 times less stringent than the ${{D}^{0}}\rightarrow
h^{\prime-}h^{+}e^{\pm}\mu^{\mp}$ limits reported in Ref. Lees  _et al._
(2020).

In this report we present a search for seven ${{D}^{0}}\rightarrow
X^{0}e^{\pm}\mu^{\mp}$ LFV decays, where $X^{0}$ represents a ${\pi}^{0}$,
${K}^{0}_{\mathrm{\scriptscriptstyle S}}$, ${\kern
1.99997pt\overline{\kern-1.99997ptK}{}}{}^{*0}$, $\rho^{0}$, $\phi$, $\omega$,
or $\eta$ meson, with data recorded with the BABAR detector at the PEP-II
asymmetric-energy ${e^{+}e^{-}}$ collider operated at the SLAC National
Accelerator Laboratory. The intermediate mesons $X^{0}$ are reconstructed
through the decays ${{\pi}^{0}}\rightarrow\gamma\gamma$,
${{K}^{0}_{\mathrm{\scriptscriptstyle S}}}\rightarrow{{\pi}^{+}}{{\pi}^{-}}$,
${{\kern
1.99997pt\overline{\kern-1.99997ptK}{}}{}^{*0}}\rightarrow{{K}^{-}}{{\pi}^{+}}$,
${{\rho}^{0}}\rightarrow{{\pi}^{+}}{{\pi}^{-}}$,
$\phi\rightarrow{{K}^{+}}\kern-1.60004pt{{K}^{-}}$,
$\omega\rightarrow{{\pi}^{+}}{{\pi}^{-}}{{\pi}^{0}}$,
$\eta\rightarrow{{\pi}^{+}}{{\pi}^{-}}{{\pi}^{0}}$, and
$\eta\rightarrow\gamma\gamma$. The branching fractions for the signal modes
are measured relative to the normalization decays
${{D}^{0}}\rightarrow{{\pi}^{-}}{{\pi}^{+}}{{\pi}^{+}}{{\pi}^{-}}$ (for
$X^{0}={{K}^{0}_{\mathrm{\scriptscriptstyle S}}},{{\rho}^{0}},\omega$),
${{D}^{0}}\rightarrow{{K}^{-}}{{\pi}^{+}}{{\pi}^{+}}{{\pi}^{-}}$
($X^{0}={{\kern 1.99997pt\overline{\kern-1.99997ptK}{}}{}^{*0}}$), and
${{D}^{0}}\rightarrow{{K}^{-}}{{K}^{+}}{{\pi}^{+}}{{\pi}^{-}}$ ($X^{0}=\phi$).
For $X^{0}={{\pi}^{0}}$ or $\eta$, the normalization mode
${{D}^{0}}\rightarrow{{K}^{-}}{{\pi}^{+}}{{\pi}^{+}}{{\pi}^{-}}$ is used as it
has the smallest branching fraction uncertainty Tanabashi  _et al._ (2018) and
the largest number of reconstructed candidates of the three normalization
modes. Although decays of the type ${{D}^{0}}\rightarrow
X^{0}h^{\prime-}h^{+}$ have momentum distributions that more closely follow
those of the signal decays under study, they suffer from smaller branching
fractions, greater uncertainties on their branching fractions, and reduced
reconstruction efficiencies relative to the three chosen normalization modes.

The ${D}^{0}$ mesons are identified using the decay
${{D}^{*+}}\rightarrow{{D}^{0}}{{\pi}^{+}}$ produced in
${e^{+}e^{-}}\rightarrow{{c}{\overline{c}}}$ events. Although ${D}^{0}$ mesons
are also produced via other processes, the use of this decay chain increases
the purity of the ${D}^{0}$ samples at the cost of a smaller number of
reconstructed ${D}^{0}$ mesons.

##  II The BABAR detector and data set

The BABAR detector is described in detail in Refs. Aubert  _et al._ (2002,
2013). Charged particles are reconstructed as tracks with a five-layer silicon
vertex detector and a 40-layer drift chamber inside a $1.5\,$T solenoidal
magnet. An electromagnetic calorimeter comprised of 6580 CsI(Tl) crystals is
used to identify and measure the energies of electrons, positrons, muons, and
photons. A ring-imaging Cherenkov detector is used to identify charged hadrons
and to provide additional lepton identification information. Muons are
primarily identified with an instrumented magnetic-flux return.

The data sample corresponds to 424$\mbox{\,fb}^{-1}$ of $e^{+}e^{-}$
collisions collected at the center-of-mass (c.m.) energy of the $\mathchar
28935\relax{(4S)}$ resonance (10.58$\mathrm{\,Ge\kern-1.00006ptV}$, on peak)
and an additional 44$\mbox{\,fb}^{-1}$ of data collected 0.04
$\mathrm{\,Ge\kern-1.00006ptV}$ below the $\mathchar 28935\relax{(4S)}$
resonance (off peak) Lees  _et al._ (2013).

Monte Carlo (MC) simulation is used to investigate sources of background
contamination and evaluate selection efficiencies. Simulated events are also
used to validate the selection procedure and for studies of systematic
effects. The signal and normalization channels are simulated with the EvtGen
package Lange (2001). We generate the signal channel decays uniformly
throughout the three-body phase space, while the normalization modes include
two-body and three-body intermediate resonances, as well as nonresonant
decays. We also generate ${e^{+}e^{-}}\rightarrow{{q}{\overline{q}}}$
($q=u,d,s,c$), Bhabha and $\mu^{+}\mu^{-}$ pairs (collectively referred to as
QED events), and ${B}{\kern 1.79993pt\overline{\kern-1.79993ptB}{}}$
background, using a combination of the EvtGen, Jetset T. Sjöstrand (1994),
KK2F Ward  _et al._ (2003), AfkQed H. Czyz and J. H. Kühn (2001), and TAUOLA
Davidson  _et al._ (2012) generators, where appropriate. The background
samples are produced with an integrated luminosity approximately 6 times that
of the data. Final-state radiation is generated using PHOTOS Golonka and Was
(2006). The detector response is simulated with GEANT 4 Agostinelli  _et al._
(2003); Allison  _et al._ (2006). All simulated events are reconstructed in
the same manner as the data.

##  III Event Selection

In the following, unless otherwise noted, all observables are evaluated in the
laboratory frame. In order to optimize the event reconstruction, candidate
selection criteria, multivariate analysis training, and fit procedure, a
rectangular area in the $m({{D}^{0}})$ versus ${\Delta
m}=m({{D}^{*+}})-m({{D}^{0}})$ plane is defined, where $m({{D}^{*+}})$ and
$m({{D}^{0}})$ are the reconstructed masses of the ${D}^{*+}$ and ${D}^{0}$
candidates, respectively. This region is kept hidden (blinded) in data until
the analysis steps are finalized. The hidden region is approximately 3 times
the root mean square (RMS) width of the $\Delta m$ and $m({{D}^{0}})$
resolutions. Its $\Delta m$ region is $0.1447<{\Delta
m}<0.1462{\mathrm{\,Ge\kern-1.00006ptV\\!/}c^{2}}$ for all modes. The
$m({{D}^{0}})$ signal peak distribution is asymmetric due to bremsstrahlung
emission, with the left-side RMS width typically $1$ to
$2$${\mathrm{\,Me\kern-1.00006ptV\\!/}c^{2}}$ wider than the right-side. The
$m({{D}^{0}})$ RMS widths vary between $5$ and
$21$${\mathrm{\,Me\kern-1.00006ptV\\!/}c^{2}}$, depending on the signal mode.

Particle identification (PID) criteria are applied to all charged daughter
tracks of the intermediate meson $X^{0}$ decays. The charged pions and kaons
are identified by measurements of their energy loss in the tracking detectors,
and the number of photons and the Cherenkov angle recorded in the ring-imaging
Cherenkov detector. These measurements are combined with information from the
electromagnetic calorimeter and the muon detector to identify electrons and
muons Aubert  _et al._ (2002, 2013). Photons are detected and their energies
are measured in the electromagnetic calorimeter. For
${{D}^{0}}\rightarrow\phi\,e^{\pm}\mu^{\mp}$, the PID requirement on the kaons
from the $\phi$ meson decay is relaxed compared to the single-kaon modes. This
increases the reconstruction efficiency for this signal mode, with little
increase in backgrounds or misidentified candidates. The muon PID requirement
depends on the signal mode, with tighter requirements imposed for modes with
more charged pions in the final state. The PID efficiency depends on the track
momentum, and is in the range $0.87-0.92$ for electrons, $0.60-0.95$ for
muons, $0.86-0.98$ for pions, and $0.84-0.92$ for kaons. The misidentification
probability Bevan  _et al._ (2014), defined as the probability that particles
are identified as one flavor (e.g. muon) that are in reality of a different
flavor (i.e. not a muon), is typically less than $0.03$ for all selection
criteria, except for the pion selection criteria, where the muon
misidentification rate can be as high as $0.35$ at low momentum.

We select events that have at least five charged tracks, except for
${{D}^{0}}\rightarrow{{\pi}^{0}}e^{\pm}\mu^{\mp}$ and
${{D}^{0}}\rightarrow\eta(\rightarrow\gamma\gamma)e^{\pm}\mu^{\mp}$, which
must have at least three. Two or more of the tracks must be identified as
leptons. The separation along the beam axis between the two leptons at their
distance of closest approach to the beam line is required to be less than
0.2$\mathrm{\,cm}$. The leptons must have opposite charges, and their momenta
must be greater than 0.3${\mathrm{\,Ge\kern-1.00006ptV\\!/}c}$. Electrons and
positrons from photon conversions are rejected by removing electron-positron
pairs with an invariant mass less than
0.03${\mathrm{\,Ge\kern-1.00006ptV\\!/}c^{2}}$ and a production vertex more
than 2$\mathrm{\,cm}$ from the beam axis.

The minimum photon energy in a signal decay is required to be greater than
0.025$\mathrm{\,Ge\kern-1.00006ptV}$. For the decays
${{D}^{0}}\rightarrow{{\pi}^{0}}e^{\pm}\mu^{\mp}$ and
${{D}^{0}}\rightarrow\eta(\rightarrow\gamma\gamma)e^{\pm}\mu^{\mp}$, the
momentum of the ${\pi}^{0}$ or $\eta$ must be greater than
0.4${\mathrm{\,Ge\kern-1.00006ptV\\!/}c}$ and the energy of each photon from
the ${\pi}^{0}$ must be greater than 0.045$\mathrm{\,Ge\kern-1.00006ptV}$. The
reconstructed ${\pi}^{0}$ invariant mass for all signal decays is required to
be between 120 and 160${\mathrm{\,Me\kern-1.00006ptV\\!/}c^{2}}$.

The reconstructed invariant masses of the ${\pi}^{0}$,
${K}^{0}_{\mathrm{\scriptscriptstyle S}}$, ${\kern
1.99997pt\overline{\kern-1.99997ptK}{}}{}^{*0}$, ${{\rho}^{0}}$, $\phi$, and
$\omega$ candidates are required to be within 19, 9, 76, 240, 20, and
34${\mathrm{\,Me\kern-1.00006ptV\\!/}c^{2}}$, of their nominal mass Tanabashi
_et al._ (2018), respectively. For the decays $\eta\rightarrow\gamma\gamma$
and $\eta\rightarrow{{\pi}^{+}}{{\pi}^{-}}{{\pi}^{0}}$, the invariant mass of
the $\eta$ candidates must be within 47 and
35${\mathrm{\,Me\kern-1.00006ptV\\!/}c^{2}}$ of the $\eta$ nominal mass,
respectively. These ranges are equivalent to 3 times the reconstructed RMS
widths.

Candidate ${D}^{0}$ mesons for the signal modes are formed from the electron
or positron, muon or antimuon, and intermediate resonance candidates. For the
normalization modes, the ${D}^{0}$ candidate is formed from four charged
tracks. Particle identification is applied to all charged tracks and the
${D}^{0}$ candidates are reconstructed with the appropriate charged-track mass
hypotheses for both the signal and normalization decays. The tracks are
required to form a good-quality vertex with a $\chi^{2}$ probability for the
vertex fit greater than 0.005. For the decay
${{D}^{0}}\rightarrow{{K}^{0}_{\mathrm{\scriptscriptstyle
S}}}e^{\pm}\mu^{\mp}$, the ${K}^{0}_{\mathrm{\scriptscriptstyle S}}$ must have
a transverse flight distance from the ${D}^{0}$ decay vertex greater than
0.2$\mathrm{\,cm}$. A bremsstrahlung energy recovery algorithm is applied to
electrons and positrons, in which the energy of photon showers that are within
a small angle (35 mrad in polar angle and 50 mrad in azimuth Aubert  _et al._
(2002)) with respect to the tangent of the initial electron or positron
direction is added to the energy of the electron or positron candidate. For
the normalization modes, the reconstructed ${D}^{0}$ meson mass is required to
be in the range
$1.81<m({{D}^{0}})<1.91{\mathrm{\,Ge\kern-1.00006ptV\\!/}c^{2}}$, while for
the signal modes, $m({{D}^{0}})$ must be in the hidden $m({{D}^{0}})$ range
defined above.

The candidate ${D}^{*+}$ is formed by combining the ${D}^{0}$ candidate with a
charged pion having a momentum greater than
$0.1{\mathrm{\,Ge\kern-1.00006ptV\\!/}c}$. For the normalization mode
${{D}^{0}}\rightarrow{{K}^{-}}{{\pi}^{+}}{{\pi}^{+}}{{\pi}^{-}}$, this pion is
required to have a charge opposite that of the kaon. The pion and ${D}^{0}$
candidate are subject to a vertex fit, with the ${D}^{0}$ mass constrained to
its known value Tanabashi  _et al._ (2018) and the requirement that the
${D}^{0}$ meson and the pion originate from the beam spot Hulsbergen (2005).
The $\chi^{2}$ probability of the fit is required to be greater than 0.005.
After the application of the ${D}^{*+}$ vertex fit, the ${D}^{0}$ candidate
momentum in the c.m. system $p^{\ast}({{D}^{0}})$ must be greater than
2.4${\mathrm{\,Ge\kern-1.00006ptV\\!/}c}$. For the normalization modes, the
mass difference $\Delta m$ is required to be $0.143<{\Delta
m}<0.148{\mathrm{\,Ge\kern-1.00006ptV\\!/}c^{2}}$, while for the signal modes
the range is $0.1395<{\Delta
m}<0.1610{\mathrm{\,Ge\kern-1.00006ptV\\!/}c^{2}}$. The extended $\Delta m$
range for the signal modes provides greater stability when fitting the
background distributions.

The requirement on the number of charged tracks strongly suppresses
backgrounds from QED processes. The $p^{\ast}({{D}^{0}})$ criterion removes
most sources of combinatorial background, as well as charm hadrons produced in
$B$ decays, which are kinematically limited to $p^{\ast}({{D}^{0}})\lesssim
2.2{\mathrm{\,Ge\kern-1.00006ptV\\!/}c}$ Aubert  _et al._ (2004).

Simulated samples indicate that the remaining background arises from
${e^{+}e^{-}}\rightarrow{{c}{\overline{c}}}$ events in which charged tracks
and neutral particles can either be lost or selected from elsewhere in the
event to form a ${D}^{0}$ candidate. To reject this background, a multivariate
selection based on a Boosted Decision Tree (BDT) discriminant is applied to
the signal modes Freund and Schapire (1997). A common set of eight input
observables is used for all modes: the momenta of the electron or positron,
muon or antimuon, and reconstructed intermediate meson; the momentum of the
lowest-momentum charged track or photon from the $X^{0}$ candidate; the
maximum angle between the direction of ${D}^{0}$ daughters and the ${D}^{0}$
direction; the total energy of all charged tracks and photons in the event,
normalized to the beam energy; the ratio
$x_{p}=p^{*}({{D}^{*+}})/\sqrt{E_{{e^{+}e^{-}}}^{*2}-m^{2}({{D}^{*+}})}$,
where $p^{*}({{D}^{*+}})$ is the c.m. momentum of the ${D}^{*+}$ candidate and
$E^{*}_{{e^{+}e^{-}}}$ is the c.m. beam energy; and the reconstructed mass of
the intermediate meson. Three additional input observables are used for the
${D}^{0}$ decays with $\omega$ or $\eta$ decaying to
${{\pi}^{+}}{{\pi}^{-}}{{\pi}^{0}}$: the momentum and reconstructed mass of
the ${\pi}^{0}$ candidate, and the energy of the lowest-energy photon from the
${\pi}^{0}$. The discriminant is trained and tested independently for each
signal mode, using simulated samples for the signal modes, and ensembles of
data outside the hidden region and
${e^{+}e^{-}}\rightarrow{{c}{\overline{c}}}$ simulated samples for the
background. Depending on the signal mode, the requirement on the discriminant
output accepts between 70% to 90% of the simulated signal sample while
rejecting between 50% to 90% of the background.

The cross feed to one signal mode from any other signal modes is estimated
from simulated samples to be less than $4\%$ in all cases, and typically less
than $1\%$, assuming equal branching fractions for all signal modes. The cross
feed to a specific normalization mode from the other two normalization modes
is predicted from simulation to be less than $0.7\%$, where the branching
fractions are taken from Ref. Tanabashi  _et al._ (2018). In the data, no
events with reconstructed normalization decays contain reconstructed signal
decays.

From the data, we find that the fraction of normalization mode events with
more than one candidate is $2.4$%, $3.6$%, and $4.4$% for
${{D}^{0}}\rightarrow{{K}^{-}}{{K}^{+}}{{\pi}^{+}}{{\pi}^{-}}$,
${{D}^{0}}\rightarrow{{K}^{-}}{{\pi}^{+}}{{\pi}^{+}}{{\pi}^{-}}$, and
${{D}^{0}}\rightarrow{{\pi}^{-}}{{\pi}^{+}}{{\pi}^{+}}{{\pi}^{-}}$,
respectively. For the signal mode with
$\eta\rightarrow{{\pi}^{+}}{{\pi}^{-}}{{\pi}^{0}}$, 40% of events have
multiple candidates. For $\eta\rightarrow\gamma\gamma$ and $\omega$ decays,
the number of events with multiple candidates is $\sim 10\%$, and for the
remaining modes it is 1% to 5%. If two or more candidates are found in an
event, the one with the highest ${D}^{*+}$ vertex $\chi^{2}$ probability is
selected. After applying the best-candidate selection, the correct ${D}^{*+}$
candidate in the simulated samples is selected with a probability of 95% or
more for the normalization modes. For the signal modes, 70% of ${D}^{*+}$
candidates are correctly selected for
$\eta\rightarrow{{\pi}^{+}}{{\pi}^{-}}{{\pi}^{0}}$, and between 86% and 94%
for the remaining modes. After the application of all selection criteria and
corrections for small differences between data and MC simulation in tracking
and PID performance, the reconstruction efficiency $\epsilon_{\rm sig}$ for
the simulated signal decays is between 1.6% and 3.6%, depending on the mode.
For the normalization decays, the reconstruction efficiency $\epsilon_{\rm
norm}$ is between $19.2$% and $24.7$%. The difference between $\epsilon_{\rm
sig}$ and $\epsilon_{\rm norm}$ is mainly due to the minimum momentum
criterion on the leptons required by the PID algorithms Aubert  _et al._
(2013).

##  IV Signal Yield Extraction

The ${{D}^{0}}\rightarrow X^{0}e^{\pm}\mu^{\mp}$ signal mode branching
fraction ${\mathcal{B}}_{\rm sig}$ is determined relative to that of the
normalization decay using

${\mathcal{B}}_{\rm sig}=\frac{N_{\rm sig}}{N_{\rm norm}}\frac{\epsilon_{\rm
norm}}{\epsilon_{\rm sig}}\frac{\mathcal{L}_{\rm norm}}{\mathcal{L}_{\rm
sig}}\frac{{\mathcal{B}}_{\rm norm}}{{\mathcal{B}}(X^{0})},$ (1)

where ${\mathcal{B}}_{\rm norm}$ is the branching fraction of the
normalization mode Tanabashi  _et al._ (2018), and $N_{\rm sig}$ and $N_{\rm
norm}$ are the fitted yields of the signal and normalization mode decays,
respectively. ${\mathcal{B}}(X^{0})$ is the branching fraction of the
intermediate meson decay channel. The symbols $\mathcal{L}_{\rm sig}$ and
$\mathcal{L}_{\rm norm}$ represent the integrated luminosities of the data
samples used for the signal ($468.2\pm 2.0$$\mbox{\,fb}^{-1}$) and the
normalization decays ($39.3\pm 0.2$$\mbox{\,fb}^{-1}$), respectively Lees  _et
al._ (2013). For the signal modes, we use both the on-peak and off-peak data
samples. For the normalization modes, a subset of the off-peak data is
sufficient for achieving statistical uncertainties that are much smaller than
the systematic uncertainties.

We perform an extended unbinned maximum likelihood fit to extract the signal
and background yields for both the normalization and signal modes Lees  _et
al._ (2014). The likelihood function is

${\mathcal{L}}=\frac{1}{N!}\exp{\left(-\sum_{j=1}^{2}n_{j}\right)}\prod_{i=1}^{N}\left[\sum_{j=1}^{2}n_{j}{\mathcal{P}}_{j}(\vec{x}_{i};\vec{\alpha}_{j})\right]\\!.$
(2)

We define the likelihood for each event candidate $i$ to be the sum of
$n_{j}{\cal P}_{j}(\vec{x}_{i};\vec{\alpha}_{j})$ over two hypotheses $j$
(signal or normalization and background). The symbol ${\cal
P}_{j}(\vec{x}_{i};\vec{\alpha}_{j})$ is the product of the probability
density functions (PDFs) for hypothesis $j$ evaluated for the measured
variables $\vec{x}_{i}$ of the $i$-th event. The total number of events in the
sample is $N$, and $n_{j}$ is the yield for hypothesis $j$. The quantities
$\vec{\alpha}_{j}$ represent parameters of ${\cal P}_{j}$. The distributions
of each discriminating variable $x_{i}$ in the likelihood function is modeled
with one or more PDFs, where the parameters $\vec{\alpha}_{j}$ are determined
from fits to signal simulation or data samples.

Each normalization mode yield $N_{\rm norm}$ is extracted by performing a two-
dimensional unbinned maximum likelihood fit to the $\Delta m$ versus
$m({{D}^{0}})$ distributions in the range $0.143<{\Delta
m}<0.148{\mathrm{\,Ge\kern-1.00006ptV\\!/}c^{2}}$ and
$1.81<m({{D}^{0}})<1.91{\mathrm{\,Ge\kern-1.00006ptV\\!/}c^{2}}$. Considering
normalization and background events separately, the measured $\Delta m$ and
$m({{D}^{0}})$ values are essentially uncorrelated and are therefore treated
as independent observables in the fits. The PDFs in the fits depend on the
normalization mode and use sums of multiple Cruijff Lees  _et al._ (2019) and
Crystal Ball Skwarnicki  functions in both $\Delta m$ and $m({{D}^{0}})$. The
functions for each observable use a common mean. The background is modeled
with an ARGUS threshold function Albrecht  _et al._ (1990) for $\Delta m$ and
a Chebyshev polynomial for $m({{D}^{0}})$. The ARGUS end point parameter is
fixed at $0.1395{\mathrm{\,Ge\kern-1.00006ptV\\!/}c^{2}}$, the $\Delta m$
kinematic threshold for ${{D}^{*+}}\rightarrow{{D}^{0}}{{\pi}^{+}}$ decays.
All other PDF parameters, together with the normalization mode and background
yields, are allowed to vary in the fit.

The fitted yields and reconstruction efficiencies for the normalization modes
are given in Table 1. Figure 1 shows projections of the unbinned maximum-
likelihood fits onto the final candidate distributions as a function of
$\Delta m$ for the normalization modes in the range $0.143<{\Delta
m}<0.148{\mathrm{\,Ge\kern-1.00006ptV\\!/}c^{2}}$.

Table 1: Summary of fitted candidate yields, with statistical uncertainties, and reconstruction efficiencies for the three normalization modes. Decay mode |  $N_{\rm norm}$ (candidates) |  $\epsilon_{\rm norm}$ (%)
---|---|---
${D}^{0}$ $\rightarrow$ ${{K}^{-}}{{\pi}^{+}}{{\pi}^{+}}{{\pi}^{-}}$ | $260\,870\pm 520$ | $20.1\pm 0.2$
${D}^{0}$ $\rightarrow$ ${{K}^{-}}{{K}^{+}}{{\pi}^{+}}{{\pi}^{-}}$ | $8480\pm 110$ | $19.2\pm 0.2$
${D}^{0}$ $\rightarrow$ ${{\pi}^{-}}{{\pi}^{+}}{{\pi}^{+}}{{\pi}^{-}}$ | $28\,470\pm 220$ | $24.7\pm 0.2$
---
Figure 1: Projections of the unbinned maximum-likelihood fits to the final
candidate distributions as a function of $\Delta m$ for the normalization
modes in the range $0.143<{\Delta
m}<0.148{\mathrm{\,Ge\kern-1.00006ptV\\!/}c^{2}}$. The solid blue line is the
total fit, the dashed red line is the signal and the dotted green line is the
background.

After the application of the selection criteria, there are on the order of 100
events or fewer available for fitting in each signal mode. Each signal mode
yield $N_{\rm sig}$ is therefore extracted by performing a one-dimensional
unbinned maximum likelihood fit to $\Delta m$ in the range $0.1395<{\Delta
m}<0.1610{\mathrm{\,Ge\kern-1.00006ptV\\!/}c^{2}}$. A Cruijff function is
implemented for the signal mode PDF, except for
${{D}^{0}}\rightarrow\phi\,e^{\pm}\mu^{\mp}$, for which two two-piece
Gaussians functions are used, and
${{D}^{0}}\rightarrow{{\rho}^{0}}e^{\pm}\mu^{\mp}$, for which two Cruijff
functions are used. The background is modeled with an ARGUS function with the
same end point used for the normalization modes. The signal PDF parameters and
the end point parameter are fixed in the fit. All other background parameters
and the signal and background yields are allowed to vary. Figure 2 shows the
results of the fits to the $\Delta m$ distributions for the signal modes.

---
Figure 2: Unbinned maximum-likelihood fits to the final candidate
distributions as a function of $\Delta m$ for the signal modes in the range
$0.1395<{\Delta m}<0.1610{\mathrm{\,Ge\kern-1.00006ptV\\!/}c^{2}}$. The solid
blue line is the total fit, the dashed red line is the signal and the dotted
green line is the background.

We test the performance of the maximum likelihood fit for the normalization
modes by generating ensembles of MC samples from the normalization and
background PDF distributions. The mean numbers of normalization and background
candidates used in the ensembles are taken from the fits to the data. The
numbers of generated background and normalization mode candidates are sampled
from a Poisson distribution. All background and normalization mode PDF
parameters are allowed to vary, except for the ARGUS function end point. No
significant biases are observed in the fitted yields of the normalization
modes. The same procedure is repeated for the maximum likelihood fits to the
signal modes, with ensembles of MC samples generated from the background PDF
distributions only, assuming a signal yield of zero. The signal PDF parameters
are fixed to the values used for the fits to the data, and the signal yield is
allowed to vary. The biases in the fitted signal yields are less than $\pm
0.3$ candidates for all modes, and these are subtracted from the fitted yields
before calculating the signal branching fractions.

To confirm the normalization procedure, the signal modes in Eq. (1) are
replaced with the decay ${{D}^{0}}\rightarrow{{K}^{-}}{{\pi}^{+}}$, which has
a well-measured branching fraction Tanabashi  _et al._ (2018). The
${{D}^{0}}\rightarrow{{K}^{-}}{{\pi}^{+}}$ decays are reconstructed using the
on-peak data sample only ($424.3\pm 1.8$$\mbox{\,fb}^{-1}$). The
${{D}^{0}}\rightarrow{{K}^{-}}{{\pi}^{+}}$ decay is selected using the same
criteria as used for the
${{D}^{0}}\rightarrow{{K}^{-}}{{\pi}^{+}}{{\pi}^{+}}{{\pi}^{-}}$ mode, which
is used as the normalization mode for this test. The
${{D}^{0}}\rightarrow{{K}^{-}}{{\pi}^{+}}$ signal yield is $1\,881\,950\pm
1380$ with $\epsilon_{\rm sig}=(27.4\pm 0.2)\%$. Thus, we determine
${\mathcal{B}}({{D}^{0}}\rightarrow{{K}^{-}}{{\pi}^{+}})=(3.98\pm 0.08\pm
0.10)$%, where the uncertainties are statistical and systematic, respectively.
This is consistent with the current world average of $(3.95\pm 0.03)$%
Tanabashi  _et al._ (2018). When the test is repeated using either
${{D}^{0}}\rightarrow{{K}^{-}}{{K}^{+}}{{\pi}^{+}}{{\pi}^{-}}$ or
${{D}^{0}}\rightarrow{{\pi}^{-}}{{\pi}^{+}}{{\pi}^{+}}{{\pi}^{-}}$ as the
normalization mode, ${\mathcal{B}}({{D}^{0}}\rightarrow{{K}^{-}}{{\pi}^{+}})$
is determined to be $(3.51\pm 0.18\pm 0.18)\%$ and $(4.12\pm 0.13\pm 0.16)\%$,
respectively.

##  V Systematic Uncertainties

The systematic uncertainties in the branching fraction determinations of the
signal modes arise from so-called additive systematic uncertainties that
affect the significance of the signal mode yields in the fits to the data
samples and from multiplicative systematic uncertainties on the luminosity and
signal reconstruction efficiencies.

The main sources of the additive systematic uncertainties in the signal yields
are associated with the model parametrizations used in the fits to the signal
modes, the fit biases, the allowed invariant-mass ranges for the ${D}^{0}$ and
$X^{0}$ candidates, the amount of cross feed, and the limited MC and data
sample sizes available for the optimization of the BDT discriminants.

The uncertainties associated with the fit model parametrizations of the signal
modes are estimated by repeating the fits with alternative PDFs. This involves
replacing the Cruijff functions with Crystal Ball functions, using a two-piece
Gaussian function, and changing the number of functions used in the PDFs. For
the background, the ARGUS function is replaced by a first- or second-order
polynomial. The largest deviation occurs when using the Crystal Ball functions
for the signal and the first-order polynomial for the background. The
systematic uncertainty is taken as half this maximum deviation. The largest
contribution comes from the normalization mode
${{D}^{0}}\rightarrow{{\pi}^{-}}{{\pi}^{+}}{{\pi}^{+}}{{\pi}^{-}}$ due to the
presence of increased background and greater uncertainty in the background
shape. To account for potential inaccuracies in the simulation of the
${D}^{0}$ and $X^{0}$ invariant mass distributions, we change the mass
selection ranges by $\pm 0.5\sigma$, where $\sigma$ is the RMS width of the
${D}^{0}$ or $X^{0}$ meson.

The systematic uncertainties in the correction on the fit biases for the
signal yields are taken from the ensembles of fits to the MC samples. Given
the central value of the signal yield obtained from the fit in each mode, the
cross feed yields from all other modes are calculated and are taken as a
systematic uncertainty. To evaluate the systematic uncertainty in the
application of the BDT discriminant, we vary the value of the selection
criterion for the BDT discriminant output, change the size of the hidden
region in data, and also retrain the BDT discriminant using a training sample
with a different ensemble of MC samples. Summing the uncertainties in
quadrature, the total additive systematic uncertainties in the signal yields
are between 0.4 and 0.9 events.

Multiplicative systematic uncertainties are due to assumptions made about the
distributions of the final-state particles in the signal simulation modeling,
the model parametrizations used in the fits to the normalization modes, the
normalization mode branching fractions, tracking and PID efficiencies, limited
simulation sample sizes, and luminosity.

Since the decay mechanism of the signal modes is unknown, we vary the angular
distributions of the simulated final-state particles from the ${D}^{0}$ signal
decay in three angular variables, defined following the prescription of Ref.
Aubert  _et al._ (2005). We weight the events, which are simulated uniformly
in phase space, using combinations of $\sin$, $\cos$, $\sin^{2}$, and
$\cos^{2}$ functions of the angular variables. The reconstruction efficiencies
calculated from simulation samples as functions of the three angles are
constant, within the statistics available. The deviations of the reweighted
efficiencies from the default average reconstruction efficiencies are
therefore small. Half the maximum change in the average reconstruction
efficiency is assigned as a systematic uncertainty.

Uncertainties associated with the fit model parametrizations of the
normalization modes are estimated by repeating the fits with alternative PDFs.
This involves swapping the Cruijff and Crystal Ball functions used in both
$\Delta m$ and $m({{D}^{0}})$. For the background, the order of the
polynomials is changed and the ARGUS function is replaced by a second-order
polynomial. Half the maximum change in the fitted yield is assigned as a
systematic uncertainty. The normalization modes branching fraction
uncertainties are taken from Ref. Tanabashi  _et al._ (2018).

For both signal and normalization modes, we include uncertainties to account
for discrepancies between reconstruction efficiencies calculated from
simulation and data samples of 1.0% per ${K}^{0}_{\mathrm{\scriptscriptstyle
S}}$, 0.8% per lepton, and 0.7% per hadron track Allmendinger  _et al._
(2013). We include a momentum-dependent ${\pi}^{0}$ reconstruction efficiency
uncertainty of 2.1% for ${{D}^{0}}\rightarrow{{\pi}^{0}}e^{\pm}\mu^{\mp}$ and
2.3% for ${{D}^{0}}\rightarrow\omega e^{\pm}\mu^{\mp}$ and
${{D}^{0}}\rightarrow\eta(\rightarrow{{\pi}^{+}}{{\pi}^{-}}{{\pi}^{0}})e^{\pm}\mu^{\mp}$.
For the PID efficiencies, we assign an uncertainty of 0.7% per track for
electrons, 1.0% for muons, 0.2% for charged pions, and 1.1% for kaons Aubert
_et al._ (2013). A systematic uncertainty of 0.4% is associated with our
knowledge of the luminosities $\mathcal{L}_{\rm norm}$ and $\mathcal{L}_{\rm
sig}$ Lees  _et al._ (2013). We assign systematic uncertainties in the range
0.8% to 1.8% to account for the limited size of the simulation samples
available for calculating reconstruction efficiencies for the signal and
normalization modes.

The simulation samples for the normalization modes contain a resonant
structure of intermediate resonances that decay to two- or three-body final
states, as well as four-body nonresonant decays. To investigate how changes in
the resonant structure affect the reconstruction efficiencies, the simulation
samples were generated using a four-body phase-space distribution only and the
reconstruction efficiencies recalculated. The resulting changes in
reconstruction efficiencies are less than the statistical uncertainties on
$\epsilon_{\rm norm}$ due to the limited size of the simulation samples, and
no systematic uncertainties are assigned. The total multiplicative systematic
uncertainties are between 4.7% and 6.8% for the normalization modes and
between 4.2% and 7.8% for the signal modes.

Table 2 summarizes the contributions of the systematic uncertainties of the
normalization modes to the systematic uncertainties in the signal mode
branching fractions, as defined in Eq. (1). Table 3 summarizes the systematic
uncertainties in the signal mode yields, excluding those due to the
normalization modes.

Table 2: Summary of the contributions to the systematic uncertainties on the signal mode branching fractions, as defined in Eq. (1), that arise from uncertainties in the measurement of the normalization modes. | ${{\pi}^{-}}{{\pi}^{+}}{{\pi}^{+}}{{\pi}^{-}}$ | ${{K}^{-}}{{\pi}^{+}}{{\pi}^{+}}{{\pi}^{-}}$ | ${{K}^{-}}{{K}^{+}}{{\pi}^{+}}{{\pi}^{-}}$
---|---|---|---
PDF variation | 4.6% | 1.0% | 1.0%
${K}^{0}_{\mathrm{\scriptscriptstyle S}}$ correction | 1.0% | 1.0% | 1.0%
Tracking correction | 3.5% | 3.5% | 3.5%
PID correction | 0.8% | 1.7% | 2.6%
Luminosity | 0.4% | 0.4% | 0.4%
Normalization $\mathcal{B}$ | 3.0% | 1.8% | 4.5%
Simulation size | 1.0% | 1.0% | 0.8%
Total |  $6.8$% |  $4.7$% |  $6.6$%
Table 3: Summary of ${{D}^{0}}\rightarrow X^{0}e^{\pm}\mu^{\mp}$ additive and multiplicative systematic uncertainties, excluding those due to the normalization modes given in Table 2. $X^{0}=$ | ${\pi}^{0}$ | ${K}^{0}_{\mathrm{\scriptscriptstyle S}}$ | ${\kern 1.99997pt\overline{\kern-1.99997ptK}{}}{}^{*0}$ | ${\rho}^{0}$ | $\phi$ | $\omega$ | $\eta$ | $\eta$
---|---|---|---|---|---|---|---|---
$X^{0}\rightarrow$ | $\gamma\gamma$ | ${{\pi}^{+}}{{\pi}^{-}}$ | ${{K}^{-}}{{\pi}^{+}}$ | ${{\pi}^{+}}{{\pi}^{-}}$ | ${{K}^{+}}{{K}^{-}}$ | ${{\pi}^{+}}{{\pi}^{-}}{{\pi}^{0}}$ | $\gamma\gamma$ | ${{\pi}^{+}}{{\pi}^{-}}{{\pi}^{0}}$
Additive (events): |  |  |  |  |  |  |
PDF variation | 0.23 | 0.05 | 0.20 | 0.16 | 0.17 | 0.26 | 0.43 | 0.16
Fit bias | 0.09 | 0.28 | 0.21 | 0.15 | 0.24 | 0.09 | 0.08 | 0.07
${D}^{0}$/$X^{0}$ mass | 0.30 | 0.04 | 0.05 | 0.07 | 0.07 | 0.07 | 0.04 | 0.23
BDT discriminant | 0.83 | 0.68 | 0.71 | 0.30 | 0.06 | 0.35 | 0.27 | 0.58
Cross feed |  |  |  | 0.01 |  |  | 0.06 |
Subtotal (candidates) | 0.92 | 0.74 | 0.76 | 0.38 | 0.31 | 0.45 | 0.52 | 0.65
Multiplicative (%): |  |  |  |  |  |  |
Angular variation | 1.4 | 2.8 | 2.0 | 3.4 | 5.3 | 1.9 | 1.6 | 1.6
${\mathcal{B}}(X^{0})$ subdecay |  | 0.1 |  |  | 1.0 | 0.8 | 0.5 | 1.2
${K}^{0}_{\mathrm{\scriptscriptstyle S}}$ correction |  | 1.0 |  |  |  |  |  |
Tracking correction | 2.3 | 3.7 | 3.7 | 3.7 | 3.7 | 3.7 | 2.3 | 3.7
PID correction | 2.7 | 2.1 | 3.0 | 2.1 | 3.9 | 3.1 | 2.7 | 3.1
${\pi}^{0}$ correction | 2.1 |  |  |  |  | 2.3 |  | 2.3
Luminosity | 0.4 | 0.4 | 0.4 | 0.4 | 0.4 | 0.4 | 0.4 | 0.4
Simulation sample size | 1.4 | 1.3 | 1.5 | 1.3 | 1.4 | 1.8 | 1.3 | 1.5
Subtotal (%) | 4.2 | 5.4 | 5.4 | 5.6 | 7.8 | 5.7 | 4.2 | 5.6

##  VI Results

Table 4 gives the fitted signal yields, reconstruction efficiencies, branching
fractions with statistical and systematic uncertainties, 90% C.L. upper limits
on the branching fractions, and previous upper limits Freyberger  _et al._
(1996); Aitala  _et al._ (2001); Tanabashi  _et al._ (2018) for the signal
modes. The yields for all the signal modes are compatible with zero. We assume
that there are no cancellations due to correlations in the systematic
uncertainties in the numerator and denominator of Eq. (1). We use the
frequentist approach of Feldman and Cousins Feldman and Cousins (1998) to
determine 90% C.L. bands. When computing the limits, the systematic
uncertainties are combined in quadrature with the statistical uncertainties in
the fitted signal yields.

Table 4: Summary of fitted signal yields $N_{\rm sig}$ with statistical and systematic uncertainties, reconstruction efficiencies $\epsilon_{\rm sig}$, branching fractions with statistical and systematic uncertainties, 90% C.L. upper limits (U.L.) on the branching fractions, and previous limits Freyberger  _et al._ (1996); Aitala  _et al._ (2001); Tanabashi  _et al._ (2018). The additive and multiplicative uncertainties are combined to obtain the overall systematic uncertainties. The branching fraction systematic uncertainties include the uncertainties in the normalization mode branching fractions. | $N_{\rm sig}$ |  | $\epsilon_{\rm sig}$ |  $\mathcal{B}$ $(\times 10^{-7})$ |  $\mathcal{B}$ 90% U.L. $(\times 10^{-7})$
---|---|---|---|---|---
Decay mode | (candidates) |  | (%) |  | BABAR | Previous
${{D}^{0}}\rightarrow{{\pi}^{0}}e^{\pm}\mu^{\mp}$ | $-0.3\pm 2.0\pm 0.9$ |  | $2.15\pm 0.03$ | $-0.6\pm 4.8\pm 2.3$ | $8.0$ | 860
${{D}^{0}}\rightarrow{{K}^{0}_{\mathrm{\scriptscriptstyle S}}}e^{\pm}\mu^{\mp}$ | $0.7\pm 1.7\pm 0.7$ |  | $3.01\pm 0.04$ | $1.9\pm 4.6\pm 1.9$ | $8.6$ | 500
${{D}^{0}}\rightarrow{{\kern 1.99997pt\overline{\kern-1.99997ptK}{}}{}^{*0}}e^{\pm}\mu^{\mp}$ | $0.8\pm 1.8\pm 0.8$ |  | $2.31\pm 0.03$ | $2.8\pm 6.1\pm 2.6$ | $12.4$ | 830
${{D}^{0}}\rightarrow{{\rho}^{0}}e^{\pm}\mu^{\mp}$ | $-0.7\pm 1.7\pm 0.4$ |  | $2.10\pm 0.03$ | $-1.8\pm 4.4\pm 1.0$ | $5.0$ | 490
${{D}^{0}}\rightarrow\phi\,e^{\pm}\mu^{\mp}$ | $0.0\pm 1.4\pm 0.3$ |  | $3.43\pm 0.04$ | $0.1\pm 3.8\pm 0.9$ | $5.1$ | 340
${{D}^{0}}\rightarrow\omega e^{\pm}\mu^{\mp}$ | $0.4\pm 2.3\pm 0.5$ |  | $1.46\pm 0.03$ | $1.8\pm 9.5\pm 1.9$ | $17.1$ | 1200
${{D}^{0}}\rightarrow\eta e^{\pm}\mu^{\mp}$ |  |  |  | $6.1\pm 9.7\pm 2.3$ | $22.5$ | 1000
with $\eta\rightarrow\gamma\gamma$ | $1.6\pm 2.3\pm 0.5$ |  | $2.96\pm 0.04$ | $7.0\pm 10.5\pm 2.4$ | $24.0$ |
with $\eta\rightarrow{{\pi}^{+}}{{\pi}^{-}}{{\pi}^{0}}$ | $0.0\pm 2.8\pm 0.7$ |  | $2.46\pm 0.04$ | $0.4\pm 25.8\pm 6.0$ | $42.8$ |

In summary, we report 90% C.L. upper limits on the branching fractions for
seven lepton-flavor-violating ${{D}^{0}}\rightarrow X^{0}e^{\pm}\mu^{\mp}$
decays. The analysis is based on a sample of ${e^{+}e^{-}}$ annihilation data
collected with the BABAR detector, corresponding to an integrated luminosity
of $468.2\pm 2.0$$\mbox{\,fb}^{-1}$. The limits are in the range
$(5.0-22.5)\times 10^{-7}$ and are between 1 and 2 orders of magnitude more
stringent than previous ${{D}^{0}}\rightarrow X^{0}e^{\pm}\mu^{\mp}$ decay
results. For the four ${{D}^{0}}\rightarrow X^{0}e^{\pm}\mu^{\mp}$ decays with
the same final state as the ${{D}^{0}}\rightarrow
h^{\prime-}h^{+}e^{\pm}\mu^{\mp}$ decays reported in Ref. Lees  _et al._
(2020), the limits are 1.5 to 3 times more stringent.

##  VII Acknowledgments

We are grateful for the extraordinary contributions of our PEP-II colleagues
in achieving the excellent luminosity and machine conditions that have made
this work possible. The success of this project also relies critically on the
expertise and dedication of the computing organizations that support BABAR.
The collaborating institutions wish to thank SLAC for its support and the kind
hospitality extended to them. This work is supported by the US Department of
Energy and National Science Foundation, the Natural Sciences and Engineering
Research Council (Canada), the Commissariat à l’Energie Atomique and Institut
National de Physique Nucléaire et de Physique des Particules (France), the
Bundesministerium für Bildung und Forschung and Deutsche
Forschungsgemeinschaft (Germany), the Istituto Nazionale di Fisica Nucleare
(Italy), the Foundation for Fundamental Research on Matter (The Netherlands),
the Research Council of Norway, the Ministry of Education and Science of the
Russian Federation, Ministerio de Economía y Competitividad (Spain), the
Science and Technology Facilities Council (United Kingdom), and the Binational
Science Foundation (U.S.-Israel). Individuals have received support from the
Marie-Curie IEF program (European Union) and the A. P. Sloan Foundation (USA).

## References

  * Paul  _et al._ (2011) A. Paul, I. I. Bigi, and S. Recksiegel, Phys. Rev. D 83, 114006 (2011).
  * Schwartz (1993) A. J. Schwartz, Mod. Phys. Lett. A 08, 967 (1993).
  * Guadagnoli and Lane (2015) D. Guadagnoli and K. Lane, Phys. Lett. B 751, 54 (2015).
  * de Boer and Hiller (2016) S. de Boer and G. Hiller, Phys. Rev. D 93, 074001 (2016).
  * Fajfer and Košnik (2015) S. Fajfer and N. Košnik, Eur. Phys. J. C 75, 567 (2015).
  * Atre  _et al._ (2009) A. Atre, T. Han, S. Pascoli, and B. Zhang, J. High Energy Phys. 05, 030 (2009).
  * Yuan  _et al._ (2013) H. Yuan, T. Wang, G.-L. Wang, W.-L. Ju, and J.-M. Zhang, J. High Energy Phys. 08, 066 (2013).
  * Hai-Rong  _et al._ (2015) D. Hai-Rong, F. Feng, and L. Hai-Bo, Chin. Phys. C 39, 013101 (2015).
  * Paul  _et al._ (2014) A. Paul, A. de la Puente, and I. I. Bigi, Phys. Rev. D 90, 014035 (2014).
  * Burdman  _et al._ (2002) G. Burdman, E. Golowich, J. A. Hewett, and S. Pakvasa, Phys. Rev. D 66, 014009 (2002).
  * Fajfer and Prelovšek (2006) S. Fajfer and S. Prelovšek, Phys. Rev. D 73, 054026 (2006).
  * Fajfer  _et al._ (2007) S. Fajfer, N. Košnik, and S. Prelovšek, Phys. Rev. D 76, 074010 (2007).
  * Aaij  _et al._ (2016) R. Aaij _et al._ (LHCb Collaboration), Phys. Lett. B 757, 558 (2016).
  * Aaij  _et al._ (2017) R. Aaij _et al._ (LHCb Collaboration), Phys. Rev. Lett. 119, 181805 (2017).
  * Lees  _et al._ (2019) J. P. Lees _et al._ (BABAR Collaboration), Phys. Rev. Lett. 122, 081802 (2019).
  * Cappiello  _et al._ (2013) L. Cappiello, O. Cata, and G. D’Ambrosio, J. High Energy Phys. 04, 135 (2013).
  * de Boer and Hiller (2018) S. de Boer and G. Hiller, Phys. Rev. D 98, 035041 (2018).
  * Kodama  _et al._ (1995) K. Kodama _et al._ (E653 Collaboration), Phys. Lett. B 345, 85 (1995).
  * Freyberger  _et al._ (1996) A. Freyberger _et al._ (CLEO Collaboration), Phys. Rev. Lett. 76, 3065 (1996).
  * Aitala  _et al._ (2001) E. M. Aitala _et al._ (E791 Collaboration), Phys. Rev. Lett. 86, 3969 (2001).
  * Ablikim  _et al._ (2018) M. Ablikim _et al._ (BESIII Collaboration), Phys. Rev. D 97, 072015 (2018).
  * Lees  _et al._ (2020) J. P. Lees _et al._ (BABAR Collaboration), Phys. Rev. Lett. 124, 071802 (2020).
  * Tanabashi  _et al._ (2018) M. Tanabashi _et al._ (Particle Data Group), Phys. Rev. D 98, 030001 (2018), and 2019 update.
  * Aubert  _et al._ (2002) B. Aubert _et al._ (BABAR Collaboration), Nucl. Instrum. Methods Phys. Res., Sect. A 479, 1 (2002).
  * Aubert  _et al._ (2013) B. Aubert _et al._ (BABAR Collaboration), Nucl. Instrum. Methods Phys. Res., Sect. A 729, 615 (2013).
  * Lees  _et al._ (2013) J. P. Lees _et al._ (BABAR Collaboration), Nucl. Instrum. Methods Phys. Res., Sect. A 726, 203 (2013).
  * Lange (2001) D. J. Lange, Nucl. Instrum. Methods Phys. Res., Sect. A 462, 152 (2001).
  * T. Sjöstrand (1994) T. Sjöstrand, Comput. Phys. Commun. 82, 74 (1994).
  * Ward  _et al._ (2003) B. F. L. Ward, S. Jadach, and Z. Was, Nucl. Phys. Proc. Suppl. 116, 73 (2003).
  * H. Czyz and J. H. Kühn (2001) H. Czyz and J. H. Kühn, Eur. Phys. J. C 18, 497 (2001).
  * Davidson  _et al._ (2012) N. Davidson, G. Nanava, T. Przedzinski, E. Richter-Was, and Z. Was, Comput. Phys. Commun. 183, 821 (2012).
  * Golonka and Was (2006) P. Golonka and Z. Was, Eur. Phys. J. C 45, 97 (2006).
  * Agostinelli  _et al._ (2003) S. Agostinelli _et al._ (GEANT 4 Collaboration), Nucl. Instrum. Methods Phys. Res., Sect. A 506, 250 (2003).
  * Allison  _et al._ (2006) J. Allison, K. Amako, J. Apostolakis, H. Araujo, P. Dubois, _et al._ (GEANT 4 Collaboration), IEEE Trans. Nucl. Sci. 53, 270 (2006).
  * Bevan  _et al._ (2014) A. J. Bevan _et al._ (BABAR and Belle Collaborations), Eur. Phys. J. C 74, 3026 (2014).
  * Hulsbergen (2005) W. D. Hulsbergen, Nucl. Instrum. Methods Phys. Res., Sect. A 552, 566 (2005).
  * Aubert  _et al._ (2004) B. Aubert _et al._ (BABAR Collaboration), Phys. Rev. D 69, 111104 (2004).
  * Freund and Schapire (1997) Y. Freund and R. E. Schapire, J. Comput. Syst. Sci. 55, 119 (1997).
  * Lees  _et al._ (2014) J. P. Lees _et al._ (BABAR Collaboration), Phys. Rev. D 89, 011102 (2014).
  * (40) T. Skwarnicki,  _A study of the radiative cascade transitions between the $\Upsilon^{\prime}$ and $\Upsilon$ resonances_, Ph.D. thesis, Institute of Nuclear Physics, Krakow, DESY-F31-86-02, 1986.
  * Albrecht  _et al._ (1990) H. Albrecht _et al._ (ARGUS Collaboration), Phys. Lett. B 241, 278 (1990).
  * Aubert  _et al._ (2005) B. Aubert _et al._ (BABAR Collaboration), Phys. Rev. D 71, 032005 (2005).
  * Allmendinger  _et al._ (2013) T. Allmendinger _et al._ , Nucl. Instrum. Methods Phys. Res., Sect. A 704, 44 (2013).
  * Feldman and Cousins (1998) G. J. Feldman and R. D. Cousins, Phys. Rev. D 57, 3873 (1998).