\section{Introduction}

The LSND \cite{lsnd} and MiniBooNE \cite{mb_oscnew,mb_oscnst} experiments have reported excesses of
$\nu_e$ and $\bar \nu_e$ charge-current quasielastic (CCQE) events in $\nu_\mu$ beams.
Exotic models beyond the three-neutrino paradigm that have been invoked to explain these anomalies include,
for example, 3+N neutrino oscillation models
involving three active neutrinos and N additional sterile neutrinos
\cite{sorel,karagiorgi,collin,giunti,giunti2,giunti3,kopp,kopp2,white_paper,3+2,diaz},
resonant neutrino oscillations \cite{resonant}, Lorentz violation
\cite{lorentz}, sterile neutrino decay \cite{sterile_decay}, scalar decay \cite{scalar_decay},
sterile neutrino nonstandard interactions
\cite{NSI}, and altered dispersion relations with sterile neutrinos \cite{sterile_extra}.
This paper presents improved MiniBooNE $\nu_e$ appearance results with increased statistics and with
additional studies that disfavor neutral-current (NC) $\pi^0$ and external event backgrounds.

\section{The MiniBooNE Experiment}

The MiniBooNE experiment makes use of the Booster Neutrino Beam (BNB) that
is produced by 8 GeV protons from
the Fermilab Booster interacting on a beryllium target inside a magnetic focusing horn,
followed by meson decay in a 50 m decay pipe.
In neutrino mode, the $\nu_\mu$, $\bar \nu_\mu$, $\nu_e$, and $\bar \nu_e$ flux
contributions at the detector are 93.5\%, 5.9\%, 0.5\%, and 0.1\%, respectively, while
in antineutrino mode, the flux
contributions are 15.7\%, 83.7\%, 0.2\%, and 0.4\%, respectively.
The $\nu_\mu$ and $\bar{\nu}_{\mu}$ fluxes peak at approximately 600 MeV and 400 MeV, respectively.
The MiniBooNE detector, described in detail in reference \cite{mb_detector}, consists of
a 12.2 m diameter sphere filled with 818 tonnes of pure mineral oil (CH$_{2}$) and
is located 541 m from the beryllium target. The detector is covered by 1520 8-inch
photomultiplier tubes (PMTs), where 1280 PMTs are in the interior detector region and 240
PMTs are located in the optically isolated outer veto region. The PMTs detect the directed
Cherenkov light and the isotropic scintillation light produced by charged particles from neutrino
interactions in the mineral oil. Events are reconstructed \cite{mb_recon} from
the hit PMT charge and time information, and the reconstructed neutrino energy, $E_\nu^{QE}$,
is estimated from the measured energy and angle of the outgoing muon
or electron, assuming the kinematics of CCQE scattering \cite{ccqe}.
The MiniBooNE experiment has collected data from 2002-2019, based on a total of $11.27 \times 10^{20}$
protons-on-target (POT) in antineutrino mode and $18.75 \times 10^{20}$ POT
in neutrino mode. Also, a special
beam off-target run collected an additional $1.86 \times 10^{20}$ POT
in a search for sub-GeV dark matter \cite{mb_dm}.
During the 17 years of
running, the BNB and MiniBooNE detector have been stable to within
3\% in neutrino energy. Fig. \ref{Michel} shows the energy distribution of Michel electrons from
stopped muon decay for the first ($6.46 \times 10^{20}$ POT from 2002 to 2007), second ($6.38 \times 10^{20}$ POT
from 2015 to 2017),
and third running periods ($5.91 \times 10^{20}$ POT from 2017 to 2019) in neutrino mode.
By adjusting the energy calibration by 2\% for the second running period and by 3\% for the third
running period, good agreement is obtained for the Michel electron energy distribution.

\begin{figure}[tbp]
\vspace{+0.1in}
\centerline{\includegraphics[angle=0, width=9.0cm]{cMichels.pdf}}
\vspace{-0.2in}
\caption{The Michel electron energy distribution
for the first, second, and third running periods in neutrino mode.
The events are normalized to the first running period. The bottom plot shows ratios of
the second and third running periods to the first running period.}
\label{Michel}
\vspace{0.1in}
\end{figure}

\section{Data Analysis}

The data analysis is optimized to measure $\nu_e$-induced CCQE
events and reject $\nu_\mu$ induced events, and is identical to the previous
analysis \cite{mb_oscnew}. Figs. \ref{PIDem}, \ref{PIDep}, and \ref{PIDpi0} show the
${\nu}_e$ CCQE data and background for the three particle identification variables in
neutrino mode in the $200<E_\nu^{QE}<1250$~MeV energy range
for the total $18.75 \times 10^{20}$ POT data. The comparison between data and background for the
full range of particle identification variables was shown in the supplementary material of
a previous publication \cite{mb_oscnst}.
The average selection efficiency is
$\sim 20\%$ ($\sim 0.1\%$) for $\nu_e$-induced CCQE events ($\nu_\mu$-induced background events)
generated over the fiducial volume.
The fraction of CCQE events in antineutrino mode that
are from wrong-sign neutrino events was determined from the angular
distributions of muons created in CCQE interactions and
by measuring CC single $\pi^+$ events \cite{wrong_sign}.
Table \ref{signal_bkgd} shows the predicted but unconstrained $\nu_e$ and
$\bar{\nu}_e$ CCQE background events for
the neutrino energy range $200<E_\nu^{QE}<1250$~MeV for both neutrino and antineutrino modes, where there
are approximately twice as many Monte Carlo events compared to data events.
Table \ref{signal_bkgd} also shows the total constrained background, where the overall normalization
of the $\nu_e$ intrinsic background is constrained by the $\nu_\mu$ CCQE event sample.
The upper limit of 1250 MeV was chosen by the collaboration before unblinding the data in 2007, while
the lower limit of 200 MeV was chosen in 2013
\cite{mb_oscnew} because it is the lowest energy for reliably reconstructing the
Cherenkov ring of $\nu_\mu$ CCQE events with a visible energy greater than 140 MeV.
From the given detector resolution estimated from the Michel electron spectrum
(Fig. \ref{Michel}), there is a negligible amount of migration from events below 200 MeV.
Note that the original lower limit was chosen to be 300 MeV before unblinding the data in 2007.
During the unblinding procedure, the lower limit was increased to 475 MeV due to the low probability
of the two-neutrino oscillation fit and worries about the single-gamma background. However, careful
studies of the single-gamma background were performed after unblinding and convinced the
collaboration that the single-gamma background was estimated correctly within systematic uncertainties
and agreed with theoretical calculations \cite{hill_zhang}.
Finally, Table \ref{signal_bkgd} shows the expected number of events corresponding to the LSND best fit
oscillation probability of 0.26\%, assuming large $\Delta m^2$ where the oscillations are washed out.
LSND and MiniBooNE have the same average value of L/E, but MiniBooNE has a larger range of L/E. Therefore,
the appearance probabilities for LSND and MiniBooNE should not be exactly the same at lower L/E values.
Figs. \ref{numuCCE} and \ref{pi0mass} show the $\nu_\mu$ CCQE $E_\nu^{QE}$ energy distribution and the
NC $\pi^0$  mass distribution in neutrino mode for the first, second, and
third running periods. As shown in the figures, the three running periods have good agreement.

\begin{figure}[tbp]
\vspace{+0.1in}
\centerline{\includegraphics[angle=0, width=9.0cm]{PIDem.pdf}}
\vspace{-0.2in}
\caption{The MiniBooNE neutrino mode
electron-muon particle identification distributions, corresponding to the total $18.75 \times 10^{20}$ POT data
in the $200<E_\nu^{QE}<1250$ MeV energy range,
for ${\nu}_e$ CCQE data (points with statistical errors) and background (colored histogram).
The dashed histogram shows the best fit to the neutrino-mode data assuming two-neutrino oscillations.}
\label{PIDem}
\vspace{0.1in}
\end{figure}

\begin{figure}[tbp]
\vspace{+0.1in}
\centerline{\includegraphics[angle=0, width=9.0cm]{PIDep.pdf}}
\vspace{-0.2in}
\caption{The MiniBooNE neutrino mode
electron-pion particle identification distributions, corresponding to the total $18.75 \times 10^{20}$ POT data
in the $200<E_\nu^{QE}<1250$ MeV energy range,
for ${\nu}_e$ CCQE data (points with statistical errors) and background (colored histogram).
The dashed histogram shows the best fit to the neutrino-mode data assuming two-neutrino oscillations.}
\label{PIDep}
\vspace{0.1in}
\end{figure}

\begin{figure}[tbp]
\vspace{+0.1in}
\centerline{\includegraphics[angle=0, width=9.0cm]{PIDpi0.pdf}}
\vspace{-0.2in}
\caption{The MiniBooNE neutrino mode
two-ring invariant-mass particle identification distributions, corresponding to the total $18.75 \times 10^{20}$ POT data
in the $200<E_\nu^{QE}<1250$ MeV energy range,
for ${\nu}_e$ CCQE data (points with statistical errors) and background (colored histogram).
The dashed histogram shows the best fit to the neutrino-mode data assuming two-neutrino oscillations.}
\label{PIDpi0}
\vspace{0.1in}
\end{figure}

\begin{table}[t]
\vspace{-0.1in}
\caption{\label{signal_bkgd} \em The expected (unconstrained) number of events
for the $200<E_\nu^{QE}<1250$~MeV neutrino
energy range from all of the backgrounds in the $\nu_e$ and $\bar{\nu}_e$
appearance analysis before using the constraint from the CC $\nu_\mu$ events.
The ``Other'' backgrounds correspond mostly to neutrino-nucleon and neutrino-electron elastic scattering.
Also shown are the constrained background, as well as
the expected number of events corresponding to the LSND best fit
oscillation probability of 0.26\%, assuming oscillations at large $\Delta m^2$.
The table shows
the diagonal-element systematic plus statistical uncertainties, which become
substantially reduced in the oscillation fits when correlations
between energy bins and between the $\nu_e$ and $\nu_\mu$ events
are included.}
\small
\begin{ruledtabular}
\begin{tabular}{ccc}
Process&Neutrino Mode&Antineutrino Mode \\
\hline
$\nu_\mu$ \& $\bar \nu_\mu$ CCQE & 107.6 $\pm$ 28.2 & 12.9 $\pm$ 4.3 \\
NC $\pi^0$ & 732.3 $\pm$ 95.5 & 112.3 $\pm$ 11.5 \\
NC $\Delta \rightarrow N \gamma$ & 251.9  $\pm$ 35.2 & 34.7 $\pm$ 5.4 \\
External Events & 109.8 $\pm$ 15.9 & 15.3 $\pm$ 2.8 \\
Other $\nu_\mu$ \& $\bar \nu_\mu$ & 130.8 $\pm$ 33.4 & 22.3 $\pm$ 3.5 \\
\hline
$\nu_e$ \& $\bar \nu_e$ from $\mu^{\pm}$ Decay & 621.1 $\pm$ 146.3 & 91.4 $\pm$ 27.6 \\
$\nu_e$ \& $\bar \nu_e$ from $K^{\pm}$ Decay & 280.7  $\pm$ 61.2 & 51.2 $\pm$ 11.0 \\
$\nu_e$ \& $\bar \nu_e$ from $K^0_L$ Decay & 79.6 $\pm$ 29.9 & 51.4 $\pm$ 18.0 \\
Other $\nu_e$ \& $\bar \nu_e$ & 8.8 $\pm$ 4.7 & 6.7 $\pm$ 6.0 \\
\hline
Unconstrained Bkgd. & $2322.6 \pm 258.3$ & $398.2  \pm 49.7$ \\
Constrained Bkgd. & $2309.4 \pm 119.6$ &  $400.6 \pm 28.5$ \\
\hline
Total Data & 2870 &478 \\
Excess & 560.6 $\pm$ 119.6 & 77.4 $\pm$ 28.5 \\
\hline
0.26\% (LSND) $\nu_\mu \rightarrow \nu_e$ & 676.3 & 100.0 \\
\end{tabular}
\vspace{-0.2in}
\end{ruledtabular}
\normalsize
\end{table}

\bigskip
\bigskip

\begin{figure}[tbp]
\vspace{+0.1in}
\centerline{\includegraphics[angle=0, width=9.0cm]{cEccqe.pdf}}
\vspace{-0.2in}
\caption{The $\nu_\mu$ CCQE muon visible energy distribution
for the first, second, and third running periods in neutrino mode.
The events are normalized to the first running period. The bottom plot shows ratios of
the second and third running periods to the first running period.}
\label{numuCCE}
\vspace{0.1in}
\end{figure}

\begin{figure}[tbp]
\vspace{+0.1in}
\centerline{\includegraphics[angle=0, width=9.0cm]{cPiMass.pdf}}
\vspace{-0.2in}
\caption{The NC $\pi^0$ mass distribution for the first, second, and third running
periods in neutrino mode. The events are normailzed to the first running period. The bottom plot shows ratios of
the second and third running periods to the first running period.}
\label{pi0mass}
\vspace{0.1in}
\end{figure}

Systematic uncertainties are determined by considering the predicted
effects on the $\nu_\mu$, $\bar{\nu}_{\mu}$, $\nu_e$, and $\bar{\nu}_e$ CCQE rates
from variations of model parameters that
include uncertainties in the neutrino and antineutrino flux estimates,
uncertainties in neutrino cross sections, uncertainties from nuclear effects, and uncertainties in
detector modeling and reconstruction.
A covariance matrix in bins of $E^{QE}_{\nu}$ is constructed
by considering the variation from each source of systematic uncertainty on the $\nu_e$ and $\bar{\nu}_e$ CCQE signal and background, and the
$\nu_\mu$ and $\bar{\nu}_{\mu}$ CCQE prediction as a function of $E_{\nu}^{QE}$.
This matrix includes correlations between any of the $\nu_e$ and $\bar{\nu}_e$ CCQE signal and background and
$\nu_\mu$ and $\bar{\nu}_{\mu}$ CCQE samples, and is used in the $\chi^2$ calculation of the oscillation fits.

\section{Electron-neutrino Appearance Results}

Figs. \ref{evis}, \ref{costh}, and \ref{excessnat} show the visible energy, $\cos \theta$, and
$E_\nu^{QE}$ distributions for
${\nu}_e$ CCQE data and background in
neutrino mode in the $200<E_\nu^{QE}<1250$~MeV energy range
for the total $18.75 \times 10^{20}$ POT data, where $\theta$ is the angle of the reconstructed
electron relative to the incident beam direction.
Each bin of reconstructed $E_\nu^{QE}$
corresponds to a distribution of ``true'' generated neutrino energies,
which can overlap adjacent bins.
In neutrino mode, a total of 2870 data events pass
the $\nu_e$ CCQE event selection requirements with $200<E_\nu^{QE}<1250$~MeV,
compared to a background expectation of $2309.4 \pm 48.1 (stat.) \pm 109.5 (syst.)$ events.
The excess, as shown in Table \ref{excesses}, is then $560.6 \pm 119.6$ events or a $4.7 \sigma$ effect.
Fig. \ref{Excess_old_new} shows the event excesses as a function of $E_\nu^{QE}$ in
neutrino mode for the first, second, and third running periods.
Combining the MiniBooNE neutrino and antineutrino data \cite{mb_oscnst}, there are a total of 3348 events
in the $200<E_\nu^{QE}<1250$~MeV energy region,
compared to a background expectation of $2710.0 \pm 52.1(stat.) \pm 122.2(syst.)$ events.
This corresponds to a total $\nu_e$ plus $\bar \nu_e$ CCQE excess of $638.0 \pm 52.1(stat.) \pm 122.2(syst.)$ events
with respect to expectation, where the statistical uncertainty is the square root of the background estimate
and the systematic uncertainty includes correlated and uncorrelated systematic parameters.
The overall significance of the excess, $4.8 \sigma$, is limited by systematic uncertainties,
assumed to be Gaussian, as the statistical significance of
the excess is $12.2 \sigma$.
The fractional unconstrained systematic
uncertainties in the $200<E_\nu^{QE}<1250$~MeV energy range are shown in Table \ref{systs}.

\begin{table}[t]
\vspace{-0.1in}
\caption{\label{excesses} \em The number of data events, background events, and excess events in neutrino
mode for different selection crteria. The errors include both statistical and systematic uncertainties.
Also shown is the significance of each event excess. R is the radius of the reconstructed event interaction
point.
}
\small
\begin{ruledtabular}
\begin{tabular}{ccccc}
Selection&Data&Background&Excess&Significance \\
\hline
$200<E_\nu^{QE}<1250$ MeV \& $R<5m$&2870&$2309.4 \pm 119.6$&$560.6 \pm 119.6$&4.7$\sigma$ \\
$150<E_\nu^{QE}<1250$ MeV \& $R<5m$&3172&$2560.4 \pm 131.5$&$611.6 \pm 131.5$&4.7$\sigma$ \\
$200<E_\nu^{QE}<1250$ MeV \& $R<4m$&1978&$1519.4 \pm 81.9$&$458.6 \pm 81.9$&5.6$\sigma$ \\
$200<E_\nu^{QE}<1250$ MeV \& $R<3m$&864&$673.9 \pm 41.2$&$190.1 \pm 41.2$&4.6$\sigma$ \\
\hline
\end{tabular}
\vspace{-0.2in}
\end{ruledtabular}
\normalsize
\end{table}

\begin{figure}[tbp]
\vspace{+0.1in}
\centerline{\includegraphics[angle=0, width=9.0cm]{nuance_datamc_data_VisibleEnergyMeV_stage_1.pdf}}
\vspace{-0.2in}
\caption{The MiniBooNE neutrino mode
visible energy distributions, corresponding to the total $18.75 \times 10^{20}$ POT data
in the $200<E_\nu^{QE}<1250$ MeV energy range,
for ${\nu}_e$ CCQE data (points with statistical errors) and background (colored histogram).
The dashed histogram shows the best fit to the neutrino-mode data assuming two-neutrino oscillations.}
\label{evis}
\vspace{0.1in}
\end{figure}

\begin{figure}[tbp]
\vspace{+0.1in}
\centerline{\includegraphics[angle=0, width=9.0cm]{nuance_datamc_data_costheta_stage_1.pdf}}
\vspace{-0.2in}
\caption{The MiniBooNE neutrino mode
$\cos \theta$ distributions, corresponding to the total $18.75 \times 10^{20}$ POT data
in the $200<E_\nu^{QE}<1250$ MeV energy range,
for ${\nu}_e$ CCQE data (points with statistical errors) and background (colored histogram).
The dashed histogram shows
the best fit to the neutrino-mode data assuming two-neutrino oscillations.}
\label{costh}
\vspace{0.1in}
\end{figure}

\begin{figure}[tbp]
\vspace{+0.1in}
\centerline{\includegraphics[angle=0, width=9.0cm]{histNu_stacked_wbf.pdf}}
\vspace{-0.2in}
\caption{The MiniBooNE neutrino mode
$E_\nu^{QE}$ distributions, corresponding to the total $18.75 \times 10^{20}$ POT data
in the $200<E_\nu^{QE}<3000$ MeV energy range,
for ${\nu}_e$ CCQE data (points with statistical errors) and predicted backgrounds (colored histograms).
The constrained background is shown as additional points with systematic error bars. The dashed histogram
shows the best fit to the neutrino-mode data assuming two-neutrino oscillations. The last bin is for the
energy interval from 1500-3000 MeV.}
\label{excessnat}
\vspace{0.1in}
\end{figure}

\begin{figure}[tbp]
\vspace{-0.0in}
\centerline{\includegraphics[angle=0, width=9.0cm]{excessNue_CompDataSets.pdf}}
\vspace{-0.2in}
\caption{
The total event excesses
in neutrino mode as a function of $E_\nu^{QE}$ for the first, second, and third running periods.
Error bars include only statistical uncertainties.}
\label{Excess_old_new}
\vspace{0.1in}
\end{figure}

\begin{table}[t]
\vspace{-0.1in}
\caption{\label{systs} \em The fractional unconstrained systematic
uncertainties in the $200<E_\nu^{QE}<1250$~MeV energy range.
}
\small
\begin{ruledtabular}
\begin{tabular}{cc}
Systematic Uncertainty&Fraction of Event Excess \\
\hline
Cross Section&35\% \\
Optical Model&23\% \\
$\pi^+$ Production&14\% \\
Neutrino Flux&7\% \\
$K^0$ Production&4\% \\
$K^+$ Production&4\% \\
\hline
\end{tabular}
\vspace{-0.2in}
\end{ruledtabular}
\normalsize
\end{table}

\begin{figure}[tbp]
\vspace{+0.1in}
\centerline{\includegraphics[angle=0, width=9.0cm]{c_hEvisCosTheta_data_data.pdf}}
\vspace{-0.2in}
\caption{The total numbers of data events in neutrino mode as functions of reconstructed
visible energy and $\cos \theta$.
There are 20 columns of visible energy from 150 to 1250 MeV and 20 rows of $\cos \theta$ from -1
to 1.}
\label{output_D}
\vspace{0.1in}
\end{figure}

\begin{figure}[tbp]
\vspace{+0.1in}
\centerline{\includegraphics[angle=0, width=9.0cm]{c_hEvisCosTheta_data_BG.pdf}}
\vspace{-0.2in}
\caption{The total numbers of background events in neutrino mode
as functions of reconstructed visible energy and $\cos \theta$.
There are 20 columns of visible energy from 150 to 1250 MeV and 20 rows of $\cos \theta$ from -1
to 1.}
\label{output_B}
\vspace{0.1in}
\end{figure}

\begin{figure}[tbp]
\vspace{+0.1in}
\centerline{\includegraphics[angle=0, width=9.0cm]{c_hEvisCosTheta_data_Excess.pdf}}
\vspace{-0.2in}
\caption{The total numbers of excess events in neutrino mode as functions of reconstructed
visible energy and $\cos \theta$.
There are 20 columns of visible energy from 150 to 1250 MeV and 20 rows of $\cos \theta$ from -1
to 1.}
\label{output_E}
\vspace{0.1in}
\end{figure}

In order to test physics models, the numbers of data events, unconstrained background events,
and excess events in neutrino mode
with visible energy between 150 and 1250 MeV are shown
in Figs. \ref{output_D}, \ref{output_B}, and \ref{output_E} as functions of visible energy and $\cos \theta$.
In these figures, there are 20 columns of visible energy from 150 to 1250 MeV and 20 rows of $\cos \theta$ from -1
to 1. There are a total of 3182 data events, 2568.8 background events and 613.2 excess events. Fig. \ref{1D_data}
shows the $\cos \theta$ distribution of data and background events
for the 20 different energy bins, while Fig. \ref{Uz_E} shows the $\cos \theta$ distributions from 0.9 to 1
for 10 different visible energy bins. Neutrino-electron elastic
scattering events are shown as the hatched region in the ``Others'' category.

Fig. \ref{Uz0.9} shows the number of data and background
events as a function of $\cos \theta$ for $\cos \theta > 0.9$, where neutrino-electron elastic
scattering events are shown as the hatched region in the ``Others'' category and contribute to the
$\cos \theta > 0.98$ bins. The neutrino-electron elastic events constitute 53\% (89\%) of the
``Others'' category for $\cos \theta > 0.90$ ($\cos \theta > 0.99$), and the category also includes
neutrino-nucleon charged-current and neutral-current scattering events.
As shown in the figure, the excess of data events over background events
is approximately the same in each bin.

Figs. \ref{EnuQE150} and \ref{Uz150} show the $E_\nu^{QE}$ and $\cos \theta$ distributions
for the $150<E_\nu^{QE}<1250$ MeV energy range, and the total event excess as a function of $E_\nu^{QE}$ is shown in
Fig. \ref{Excess_new}.
The solid curve on the latter plot shows the two-neutrino oscillation prediction at the
best-fit point ($\sin^22\theta = 0.807$, $\Delta m^2=0.043$ eV$^2$).
The lowest energy data point has less acceptance than the other data points due to the requirement
that the visible energy be greater than 140 MeV.
Table \ref{excesses} lists the number of data events, background events, excess events, and
excess significance for the $150<E_\nu^{QE}<1250$ MeV energy range.

\begin{figure}[tbp]
\vspace{+0.1in}
\centerline{\includegraphics[angle=0, width=18.0cm]{c1DX_h1DX_hEvisCosTheta_data_0.pdf}}
\vspace{-0.2in}
\caption{The $\cos \theta$ distribution of data events (points with error bars)
and background events (colored histogram) in neutrino mode for the 20 different visible energy bins
from 150 to 1250 MeV.}
\label{1D_data}
\vspace{0.1in}
\end{figure}

\begin{figure}[tbp]
\vspace{+0.1in}
\centerline{\includegraphics[angle=0, width=18.0cm]{c1DX_h1DX_hEvisCosTheta_09_data_0.pdf}}
\vspace{-0.2in}
\caption{The $\cos \theta$ distributions from 0.9 to 1 of data events (points with statistical errors)
and background events (histogram) in neutrino mode for 10 different visible energy bins
from 150 to 1250 MeV. Neutrino-electron elastic scattering events are shown as the hatched
region in the ``Others'' category.}
\label{Uz_E}
\vspace{0.1in}
\end{figure}

\begin{figure}[tbp]
\vspace{+0.1in}
\centerline{\includegraphics[angle=0, width=9.0cm]{cCosTheta.pdf}}
\vspace{-0.2in}
\caption{The MiniBooNE neutrino mode
$\cos \theta$ distribution for $\cos \theta > 0.9$, corresponding to the total $18.75 \times 10^{20}$ POT neutrino data
in the visible energy range from 150 to 1250 MeV,
for ${\nu}_e$ CCQE data (points with statistical errors) and predicted backgrounds (colored histograms).
Neutrino-electron elastic scattering events are shown as the hatched region in the ``Others'' category.}
\label{Uz0.9}
\vspace{0.1in}
\end{figure}

\begin{figure}[tbp]
\vspace{+0.1in}
\centerline{\includegraphics[angle=0, width=9.0cm]{EnuQE_150_1250.pdf}}
\vspace{-0.2in}
\caption{The MiniBooNE neutrino mode
$E_\nu^{QE}$ distributions, corresponding to the total $18.75 \times 10^{20}$ POT neutrino data
in the $150<E_\nu^{QE}<1250$ MeV energy range,
for ${\nu}_e$ CCQE data (points with statistical errors) and predicted backgrounds (colored histograms).
The dashed histogram
shows the best fit to the neutrino-mode data assuming two-neutrino oscillations.}
\label{EnuQE150}
\vspace{0.1in}
\end{figure}

\begin{figure}[tbp]
\vspace{+0.1in}
\centerline{\includegraphics[angle=0, width=9.0cm]{Uz_150_1250.pdf}}
\vspace{-0.2in}
\caption{The MiniBooNE neutrino mode
$\cos \theta$ distributions, corresponding to the total $18.75 \times 10^{20}$ POT neutrino data
in the $150<E_\nu^{QE}<1250$ MeV energy range,
for ${\nu}_e$ CCQE data (points with statistical errors) and predicted backgrounds (colored histograms).
The dashed histogram
shows the best fit to the neutrino-mode data assuming two-neutrino oscillations.}
\label{Uz150}
\vspace{0.1in}
\end{figure}

\begin{figure}[tbp]
\vspace{-0.0in}
\centerline{\includegraphics[angle=0, width=18.0cm]{cBestFitSigma-2.pdf}}
\vspace{-0.2in}
\caption{
The total event excess
in neutrino mode, corresponding to $18.75 \times 10^{20}$ POT in the $150<E_\nu^{QE}<3000$ MeV energy range.
The solid curve shows the best fit to the neutrino-mode and antineutrino-mode data
in the $200 < E_\nu^{QE} < 1250$ MeV energy range
assuming two-neutrino oscillations, while the dashed curves show selected points on the 1 sigma contour.
Also shown is the 1-sigma allowed band.
The outer error bars include statistical plus constrained systematic uncertainties, while the inner
error bars show the statistical uncertainties. The lowest energy data point
shows only the unconstrained background statistical uncertainty.}
\label{Excess_new}
\vspace{0.1in}
\end{figure}

\section{Neutrino Oscillation Fits}

Fig.~\ref{limitab2} shows the MiniBooNE allowed regions in both neutrino mode and antineutrino
mode \cite{mb_oscnst}
for events with $200 < E^{QE}_{\nu} < 3000$ MeV within a two-neutrino oscillation model.
For this oscillation fit the entire data set is used and includes the $18.75 \times 10^{20}$ POT
data in neutrino mode and the $11.27 \times 10^{20}$ POT data in antineutrino mode.
Also shown are 90\% C.L. limits from the KARMEN \cite{karmen}
and OPERA \cite{opera} experiments.
The best combined neutrino oscillation fit occurs at
($\sin^22\theta$, $\Delta m^2$) $=$ (0.807, 0.043 eV$^2$).
The $\chi^2/ndf$ for the best-fit point in the energy range $200<E_\nu^{QE}<1250$~MeV
is 21.7/15.5 with a probability of $12.3\%$, and the background-only fit
has a  $\chi^2$ probability of $3 \times 10^{-7}$ relative to the best oscillation fit
and a $\chi^2/ndf = 50.7/17.3$ with a probability of $0.01\%$.

\begin{figure}[tbp]
\vspace{-0.25in}
 \centerline{\includegraphics[width=9.0cm]{miniboone_fake_oct19_contNunubar_with_karmen.pdf}}
 \vspace{-0.3in}
\caption{MiniBooNE allowed regions for combined neutrino mode ($18.75
  \times 10^{20}$ POT) and antineutrino mode ($11.27 \times 10^{20}$ POT)
data sets for events with
$200 < E^{QE}_{\nu} < 3000$ MeV within a two-neutrino oscillation model.
The shaded areas show the 90\% and 99\% C.L. LSND
$\bar{\nu}_{\mu}\rightarrow\bar{\nu}_e$ allowed
regions. The black point shows the MiniBooNE best fit point. Also shown are 90\% C.L. limits
from the KARMEN \cite{karmen}
and OPERA \cite{opera} experiments.}
\label{limitab2}
\vspace{0.1in}
\end{figure}

Fig. \ref{LoverE} compares the $L/E_\nu^{QE}$
distributions for the MiniBooNE data excesses in neutrino mode and
antineutrino mode to the $L/E$ distribution from LSND \cite{lsnd}.
The error bars show statistical uncertainties only. As shown in the figure,
there is agreement among all three data sets.
Assuming two-neutrino oscillations, the curves show fits to the
MiniBooNE data described above.
The significance of the combined LSND ($3.8 \sigma$) \cite{lsnd}
and MiniBooNE ($4.8 \sigma$) excesses is $6.1 \sigma$, which is obtained by adding the significances
in quadrature, as the two experiments have completely different neutrino energies, neutrino fluxes,
reconstructions, backgrounds, and systematic uncertainties.

\begin{figure}[tbp]
\vspace{-0.0in}
\centerline{\includegraphics[angle=0, width=9.0cm]{LoverE.pdf}}
\vspace{-0.2in}
\caption{A comparison between the $L/E_\nu^{QE}$
distributions for the MiniBooNE data excesses in neutrino mode
($18.75 \times 10^{20}$ POT) and antineutrino mode ($11.27
\times 10^{20}$ POT) to the $L/E$ distribution from LSND \cite{lsnd}.
The error bars show
statistical uncertainties only. The curves show fits to the
MiniBooNE data, assuming two-neutrino oscillations, while the shaded area is the MiniBooNE $1 \sigma$
allowed band.
The best-fit curve corresponds to ($\sin^22\theta$, $\Delta m^2$) $=$ (0.807, 0.043 eV$^2$), while
the dashed curve corresponds to a $1 \sigma$ fit point
at ($\sin^22\theta$, $\Delta m^2$) $=$ (0.01, 0.4 eV$^2$).}
\label{LoverE}
\vspace{0.1in}
\end{figure}

\section{Background Studies and Constraints}

Constraints have been placed on the various backgrounds in Table \ref{signal_bkgd} by direct
measurements of these backgrounds in the MiniBooNE detector.
The $\nu_\mu$ CC background has been well measured \cite{mb_numuccqe}
by using the Michel electrons from muon decay to identify
the event topology. Likewise, the NC $\pi^0$ background has also been well measured \cite{mb_pi0}
by reconstructing the
two-gamma invariant mass.

In addition, a fit to the vertex radial distribution, shown in
Fig. \ref{Rdist}, allows a
constraint to be placed on the NC $\pi^0$ background, due to this background having more events near the edge of the
5 m radius fiducial volume. (NC $\pi^0$ events near the edge of the fiducial volume have a greater chance of one
photon leaving the detector with the remaining photon then mis-reconstructing as an electron candidate.)
Fig. \ref{RExcess} shows the excess event radial distributions, where
different processes are normalized to explain the event excess, while
Table \ref{Rfit} shows the result of log-likelihood shape-only fits to the radial
distribution and the multiplicative factor that is required for each hypothesis
to explain the observed event excess. The two-neutrino oscillation hypothesis
fits the radial distribution best with a $\chi^2 = 8.4/9 ndf$, while the NC $\pi^0$ hypothesis has a worse fit
with a $\chi^2 = 17.2/9 ndf$.
The intrinsic $\nu_e$ backgrounds have a worse $\chi^2$ than the two-neutrino
oscillation hypothesis due to higher energy $\nu_e$ events having a different radial distribution than lower
energy $\nu_e$ events.

\begin{figure}[tbp]
\vspace{-0.25in}
 \centerline{\includegraphics[angle=0,width=9.0cm]{nuance_datamc_data_R500cm3_stage_1.pdf}}
 \vspace{-0.3in}
\caption{The MiniBooNE radial vertex distribution, corresponding to the total $18.75 \times 10^{20}$
POT data in neutrino mode in the $200<E_\nu^{QE}<1250$ MeV energy range,
for ${\nu}_e$ CCQE data (points with statistical errors) and background (histogram).
The dashed histogram shows the best fit to the neutrino-mode data assuming two-neutrino oscillations.}
\label{Rdist}
\vspace{0.1in}
\end{figure}

\begin{figure}[tbp]
\vspace{-0.25in}
 \centerline{\includegraphics[angle=0,width=18.0cm]{rshapes_v2.pdf}}
 \vspace{-0.3in}
\caption{The excess event radial distributions in neutrino mode with only statistical errors
in the $200<E_\nu^{QE}<1250$ MeV energy range, where
different processes are normalized to explain the event excess. The different processes are the
following: (a) $\Delta \rightarrow N \gamma$; (b) External Events; (c) $\nu_e$ \& $\bar \nu_e$ from $K^0_L$ Decay;
(d) $\nu_e$ \& $\bar \nu_e$ from $K^{\pm}$ Decay; (e) $\nu_e$ \& $\bar \nu_e$ from $\mu^{\pm}$ Decay;
(f) Other $\nu_e$ \& $\bar \nu_e$; (g) NC $\pi^0$; (h) Best Fit oscillations.}
\label{RExcess}
\vspace{0.1in}
\end{figure}

\begin{figure}[tbp]
\vspace{-0.25in}
 \centerline{\includegraphics[angle=0,width=9.0cm]{h1_TRWMZb1MeVBinHist.pdf}}
 \vspace{-0.3in}
\caption{The bunch timing for data events in neutrino mode compared to the expected background
in the $200<E_\nu^{QE}<1250$ MeV energy range. Almost all of the excess data
events occur, as expected from neutrino events in the detector, within the first 8 ns of the bunch timing.
This data sample uses events collected with the new fiber timing system and represents about 40\% of the entire neutrino mode sample.}
\label{Timing}
\vspace{0.1in}
\end{figure}

\bigskip

\begin{table}[t]
\vspace{-0.1in}
\caption{\label{Rfit} \em The result of log-likelihood shape-only fits to the radial
distribution in neutrino mode,
assuming only statistical errors, where different processes are normalized to explain the observed event excess.
The two-neutrino hypothesis
fits the radial distribution best with a $\chi^2 = 8.4/9 ndf$, while the NC $\pi^0$ hypothesis has a worse fit
with a $\chi^2 = 17.2/9 ndf$. Also shown is the multiplicative factor that is required for each hypothesis
to explain the observed event excess.
}
\small
\begin{ruledtabular}
\begin{tabular}{ccc}
Hypothesis&Multiplicative factor&$\chi^2 /9 ndf$ \\
\hline
NC $\Delta \rightarrow N \gamma$ Background & 3.18 & 10.0 \\
External Event Background & 5.98 & 44.9 \\
$\nu_e$ \& $\bar \nu_e$ from $K^0_L$ Decay Background & 7.85 & 14.8 \\
$\nu_e$ \& $\bar \nu_e$ from $K^{\pm}$ Decay Background & 2.95 & 16.3 \\
$\nu_e$ \& $\bar \nu_e$ from $\mu^{\pm}$ Decay Background & 1.88 & 16.1 \\
Other $\nu_e$ \& $\bar \nu_e$ Background & 3.21 & 12.5 \\
NC $\pi^0$ Background & 1.75 & 17.2 \\
Best Fit Oscillations & 1.24 & 8.4 \\
\hline
\end{tabular}
\vspace{-0.2in}
\end{ruledtabular}
\normalsize
\end{table}

\bigskip
\bigskip

Single-gamma backgrounds
from external neutrino interactions (``dirt" backgrounds) are estimated using topological and spatial cuts
to isolate the events whose vertices are near the edge of the detector and
point towards the detector center \cite{mb_lowe}.
The external event background estimate has been confirmed by measuring the absolute time of signal
events relative to the proton beam microstructure (52.81 MHz extraction frequency), which corresponds to buckets of beam approximately
every 18.9 ns.
Fig. \ref{Timing} shows that the event excess peaks in the 8 ns window associated with beam bunch time,
as expected from neutrino events in the detector, and is inconsistent with external neutrino events or
beam-off events, which would be approximately flat in time. Also, the observed background level outside
of the beam agrees well with the predicted background estimate. In addition, good agreement is obtained for the
event excess with $\cos \theta >0.9$.
The timing reconstruction performed here is similar to the reconstruction in
reference \cite{mb_dm}, but with a different time offset applied.

The $\Delta \rightarrow N + \gamma$ background is determined from the NC $\pi^0$ event sample \cite{mb_pi0}, which
has contributions from $\Delta$ production in $^{12}$C (52.2\%), $\Delta$ production in H$_2$ (15.1\%),
coherent scattering on $^{12}$C (12.5\%), coherent scattering on H$_2$ (3.1\%), higher-mass resonances
(12.9\%), and non-resonant background (4.2\%). The fraction of $\Delta$ decays to $\pi^0$
is 2/3 from the Clebsch-Gordon coefficients, and the probability of pion escape from the $^{12}$C nucleus
is estimated to be 62.5\%.
The $\Delta$ radiative branching fraction is 0.60\% for $^{12}$C and 0.68\% for H$_2$
after integration over all the invariant mass range,
where the single gamma production branching ratio increases below the pion production threshold.
With these values, the ratio of single gamma events to NC $\pi^0$ events, $R$, can
be estimated to be $$R=0.151 \times 0.0068 \times 1.5 + 0.522 \times 0.0060 \times 1.5 /0.625 = 0.0091.$$
Note that single gamma events are assumed to come entirely from $\Delta$ radiative decay. The total uncertainty
on this ratio is 14.0\% (15.6\%) in neutrino (antineutrino) mode. This estimate of $R=0.0091 \pm 0.0013$
agrees fairly well with theoretical calculations of the single gamma event rate \cite{hill_zhang}.

The intrinsic $\nu_e$ background comes almost entirely from muon and kaon
decay-in-flight in the beam decay pipe.
MiniBooNE $\nu_\mu$ CCQE event measurements \cite{mb_numuccqe}
constrain the size and energy dependence of the intrinsic $\nu_e$ background from muon decay,
while the intrinsic $\nu_e$ background from kaon decay is constrained by fits to kaon production
data and SciBooNE measurements \cite{sciboone_kaon}. Furthermore, due to the higher energy of the intrinsic
$\nu_e$ background, this background is disfavored from the fit
to the radial distribution, as shown in Table \ref{Rfit}.

Finally, backgrounds from exotic $\pi^0$ decay in the neutrino production target are ruled out from
the MiniBooNE beam-dump run, where the incident proton beam was steered above the Be target and
interacted in the steel beam dump at the downstream end of the decay pipe. No excess of events
was observed \cite{mb_dm}, which set limits on light dark matter and other exotic $\pi^0$ decays.

\section{Constraints on NC $\pi^0$ Background with Tighter Radius Selection}

Explanations for the event excess have included unsimulated photons entering the detector from external interactions and
the undersimulation of photons lost from $\pi^0$ production within the detector.   To test these explanations in a model-independent way,
we can use our higher event statistics to study the change in the excess as a function
of tighter fiducial volume cuts. The NC $\pi^0$ and external event backgrounds preferentially
populate higher radius compared to electron neutrino interactions. Therefore, reducing the fiducial radius is expected to
reduce the significance of the excess if it is due to these backgrounds and increase the significance of
the excess if its distribution is $\nu_e$-like. If we change the standard 5 m cut to 4 m, we find
there are 1978 data events in neutrino mode, $1519.4 \pm 81.9$ background events, and an excess of
$458.6 \pm 81.9$ events ($5.6 \sigma$). If we
use a 3 m cut, we find 864 data events, $673.9 \pm 41.2$ background events, and
an excess of $190.1 \pm 41.2$ events ($4.6 \sigma$), consistent with what is expected
if the signal is more $\nu_e$-like. The event statistics are shown in Table \ref{excesses},
while Figs. \ref{EnuQE_Rlt4} and \ref{Uz_Rlt4}
show the reconstructed neutrino energy and $\cos \theta$ distributions for electronlike events with radius less than 4 m for
both data events and background events.

\begin{figure}[tbp]
\vspace{-0.25in}
 \centerline{\includegraphics[angle=0,width=9.0cm]{nuance_datamc_data_E_QEMeV_stage_1_lt400cm.pdf}}
 \vspace{-0.3in}
\caption{The MiniBooNE neutrino mode
$E_\nu^{QE}$ distributions, corresponding to the total $18.75 \times 10^{20}$ POT data in neutrino mode
in the $200<E_\nu^{QE}<1250$ MeV energy range,
for ${\nu}_e$ CCQE data (points with statistical errors) and background
(histogram) with radius less than 4 m. The dashed histogram shows
the best fit to the neutrino-mode data assuming two-neutrino oscillations.}
\label{EnuQE_Rlt4}
\vspace{0.1in}
\end{figure}

\begin{figure}[tbp]
\vspace{-0.25in}
 \centerline{\includegraphics[angle=0,width=9.0cm]{nuance_datamc_data_costheta_stage_1_lt400cm.pdf}}
 \vspace{-0.3in}
\caption{The MiniBooNE neutrino mode
$\cos \theta$ distributions, corresponding to the total $18.75 \times 10^{20}$ POT data in neutrino mode
in the $200<E_\nu^{QE}<1250$ MeV energy range,
for ${\nu}_e$ CCQE data (points with statistical errors) and background
(histogram) with radius less than 4 m. The dashed histogram shows
the best fit to the neutrino-mode data assuming two-neutrino oscillations.}
\label{Uz_Rlt4}
\vspace{0.1in}
\end{figure}

\section{Conclusion}

In summary, the MiniBooNE experiment
observes a total excess
of $638.0 \pm 52.1(stat.) \pm 132.8(syst.)$ electronlike events in the energy range $200<E_\nu^{QE}<1250$~MeV
in both neutrino and antineutrino running modes.
The overall significance of the excess, $4.8 \sigma$, is limited by systematic uncertainties,
assumed to be Gaussian, as the statistical significance of
the excess is $12.2 \sigma$.
All of the major backgrounds are constrained by in situ
event measurements.
Beam timing information shows that almost all of the excess is in time with neutrinos that interact
in the detector. The radius distribution shows that the excess is distributed throughout the volume,
while tighter cuts on the fiducal volume increase the significance of the excess.
The data likelihood ratio disfavors models that
explain the event excess due to entering or exiting photons.
The MiniBooNE event excess will be further studied by the Fermilab
short-baseline neutrino (SBN) program \cite{sbn} and by the JSNS$^2$ experiment at J-PARC \cite{jsns2}.

\begin{acknowledgments}
We acknowledge the support of Fermilab, the Department of Energy,
and the National Science Foundation, and
we acknowledge Los Alamos National Laboratory for LDRD funding.
\vspace{-0.01in}
\end{acknowledgments}

\vspace{-0.2in}


\section*{Acknowledgements}
\noindent We express our gratitude to our colleagues in the CERN
accelerator departments for the excellent performance of the LHC. We
thank the technical and administrative staff at the LHCb
institutes.
We acknowledge support from CERN and from the national agencies:
CAPES, CNPq, FAPERJ and FINEP (Brazil);
MOST and NSFC (China);
CNRS/IN2P3 (France);
BMBF, DFG and MPG (Germany);
INFN (Italy);
NWO (Netherlands);
MNiSW and NCN (Poland);
MEN/IFA (Romania);
MSHE (Russia);
MinECo (Spain);
SNSF and SER (Switzerland);
NASU (Ukraine);
STFC (United Kingdom);
DOE NP and NSF (USA).
We acknowledge the computing resources that are provided by CERN, IN2P3
(France), KIT and DESY (Germany), INFN (Italy), SURF (Netherlands),
PIC (Spain), GridPP (United Kingdom), RRCKI and Yandex
LLC (Russia), CSCS (Switzerland), IFIN-HH (Romania), CBPF (Brazil),
PL-GRID (Poland) and OSC (USA).
We are indebted to the communities behind the multiple open-source
software packages on which we depend.
Individual groups or members have received support from
AvH Foundation (Germany);
EPLANET, Marie Sk\l{}odowska-Curie Actions and ERC (European Union);
A*MIDEX, ANR, Labex P2IO and OCEVU, and R\'{e}gion Auvergne-Rh\^{o}ne-Alpes (France);
Key Research Program of Frontier Sciences of CAS, CAS PIFI, and the Thousand Talents Program (China);
RFBR, RSF and Yandex LLC (Russia);
GVA, XuntaGal and GENCAT (Spain);
the Royal Society
and the Leverhulme Trust (United Kingdom).

\section{${\ensuremath{{\PJ\mskip -3mu/\mskip -2mu\Ppsi}}}\xspace$-pair mass distributions in bins of ${\ensuremath{p^{{\rm di}{\text -}\PJ\mskip -3mu/\mskip -2mu\Ppsi}_{\rm T}}}\xspace$}
\unboldmath
\label{app:Mbinned}
\begin{figure}[htbp]
\begin{center}
\includegraphics[width=1.0\linewidth]{./figs/FigS1.pdf}
\vspace*{-0.1cm}
\end{center}
   \caption{
   Invariant mass spectra of ${\ensuremath{{\PJ\mskip -3mu/\mskip -2mu\Ppsi}}}\xspace$-pair candidates in the six ${\ensuremath{p^{{\rm di}{\text -}\PJ\mskip -3mu/\mskip -2mu\Ppsi}_{\rm T}}}\xspace$ regions with boundaries
   $\left\{ 0, 5, 6, 8, 9.5, 12, 50 \right\}\ensuremath{\aunit{Ge\kern -0.1em V\!/}c}\xspace$
   with reconstructed ${\ensuremath{{\PJ\mskip -3mu/\mskip -2mu\Ppsi}}}\xspace$ masses in the  (black)~signal and (blue)~background regions,
   respectively.
   }\label{fig:Xmassbinned}
\end{figure}

\boldmath
\section{Additional fits to the ${\ensuremath{{\PJ\mskip -3mu/\mskip -2mu\Ppsi}}}\xspace$-pair mass spectrum}
\unboldmath
\label{app:morefits}
Figure~\ref{fig:fitXcut2} shows the fits to the ${\ensuremath{{\PJ\mskip -3mu/\mskip -2mu\Ppsi}}}\xspace$-pair mass spectrum with (a) the threshold structure described by a single Breit--Wigner (BW) lineshape and (b) using a model that contains a single BW resonance interfering with the SPS continuum.
The $\ensuremath{\chi^2}\xspace/{\rm ndof}$ of the two fits
are $125.6/92$
and $118.6/91$,
corresponding to a probability of $1.2\%$ and $2.8\%$, respectively.
Figure~\ref{fig:fitXcut3} shows the fit with an additional BW function introduced to describe the $7.2\ensuremath{\aunit{Ge\kern -0.1em V\!/}c^2}\xspace$ structure, based on the model that contains two BW lineshapes for the threshold structure and a BW shape for the $X(6900)$ structure on top of the NRSPS plus DPS continuum.

\begin{figure}[h]
\begin{center}
\sidesubfloat[]{\includegraphics[width=0.48\linewidth]{./figs/FigS2a.pdf}}
\sidesubfloat[]{\includegraphics[width=0.48\linewidth]{./figs/FigS2b.pdf}}\\
\vspace*{-0.1cm}
\end{center}
   \caption{Invariant mass spectra of weighted di-${\ensuremath{{\PJ\mskip -3mu/\mskip -2mu\Ppsi}}}\xspace$ candidates with ${\ensuremath{p^{{\rm di}{\text -}\PJ\mskip -3mu/\mskip -2mu\Ppsi}_{\rm T}}}\xspace>5.2\ensuremath{\aunit{Ge\kern -0.1em V\!/}c}\xspace$ and overlaid projections of the ${\ensuremath{p^{{\rm di}{\text -}\PJ\mskip -3mu/\mskip -2mu\Ppsi}_{\rm T}}}\xspace$-threshold fit with (a)~the threshold structure described as a single BW function,
   and (b)~assuming a single BW
   interfering with the SPS continuum.
   }\label{fig:fitXcut2}
\end{figure}

\begin{figure}[h]
\begin{center}
\includegraphics[width=0.48\linewidth]{./figs/FigS3.pdf}
\vspace*{-0.1cm}
\end{center}
   \caption{Invariant mass spectra of weighted di-${\ensuremath{{\PJ\mskip -3mu/\mskip -2mu\Ppsi}}}\xspace$ candidates with ${\ensuremath{p^{{\rm di}{\text -}\PJ\mskip -3mu/\mskip -2mu\Ppsi}_{\rm T}}}\xspace>5.2\ensuremath{\aunit{Ge\kern -0.1em V\!/}c}\xspace$ and overlaid projections of the ${\ensuremath{p^{{\rm di}{\text -}\PJ\mskip -3mu/\mskip -2mu\Ppsi}_{\rm T}}}\xspace$-threshold fit with an additional BW function introduced to describe the $7.2\ensuremath{\aunit{Ge\kern -0.1em V\!/}c^2}\xspace$ structure, based on the model that contains two BW lineshapes for the threshold structure and a BW shape for the $X(6900)$ structure on top of the NRSPS plus DPS continuum.
   }\label{fig:fitXcut3}
\end{figure}

\boldmath
\section{Supplement to fit result of model II}
\unboldmath
\label{app:itffit}
In model II that contains a BW lineshape for the threshold structure interfering with the NRSPS, a BW shape for the $X(6900)$ structure and the DPS continuum,
the parameters of the lower-mass BW lineshape is determined to $M=6741\pm6\aunit{(stat)}\xspace\ensuremath{\aunit{Me\kern -0.1em V\!/}c^2}\xspace$ and $\Gamma=288\pm16\aunit{(stat)}\xspace\aunit{Me\kern -0.1em V}\xspace$.
The systematic uncertainties on the mass and natural width are not studied.
Due to the complex nature of the threshold structure, and the simple interference scenario considered,
this study is not considered to claim a state with the parameters reported here.

Projections of the fit to the ${\ensuremath{{\PJ\mskip -3mu/\mskip -2mu\Ppsi}}}\xspace$-pair invariant mass spectra in bins of ${\ensuremath{p^{{\rm di}{\text -}\PJ\mskip -3mu/\mskip -2mu\Ppsi}_{\rm T}}}\xspace$ assuming the interference between the threshold structure and the SPS continuum are shown in Fig.~\ref{fig:fitXbinned_itf}.

\begin{figure}[!htbp]
\begin{center}
   \sidesubfloat[]{\includegraphics[width=0.45\linewidth]{./figs/FigS4a.pdf}}
   \sidesubfloat[]{\includegraphics[width=0.45\linewidth]{./figs/FigS4b.pdf}}\\
   \sidesubfloat[]{\includegraphics[width=0.45\linewidth]{./figs/FigS4c.pdf}}
   \sidesubfloat[]{\includegraphics[width=0.45\linewidth]{./figs/FigS4d.pdf}}\\
   \sidesubfloat[]{\includegraphics[width=0.45\linewidth]{./figs/FigS4e.pdf}}
   \sidesubfloat[]{\includegraphics[width=0.45\linewidth]{./figs/FigS4f.pdf}}\\
\vspace*{-0.1cm}
\end{center}
   \caption{Invariant mass spectra of weighted di-$/{\ensuremath{{\PJ\mskip -3mu/\mskip -2mu\Ppsi}}}\xspace$ candidates in bins of ${\ensuremath{p^{{\rm di}{\text -}\PJ\mskip -3mu/\mskip -2mu\Ppsi}_{\rm T}}}\xspace$ and overlaid projections of the ${\ensuremath{p^{{\rm di}{\text -}\PJ\mskip -3mu/\mskip -2mu\Ppsi}_{\rm T}}}\xspace$-binned fit assuming that the threshold structure interferes with the SPS continuum.
   }\label{fig:fitXbinned_itf}
\end{figure}

\boldmath
\section{Dependence of experimental efficiency on ${\ensuremath{M_{{\rm di}{\text -}\PJ\mskip -3mu/\mskip -2mu\Ppsi}}}\xspace$}
\unboldmath
The variation of the experimental efficiency with respect to ${\ensuremath{M_{{\rm di}{\text -}\PJ\mskip -3mu/\mskip -2mu\Ppsi}}}\xspace$ is shown in Fig.~\ref{fig:eff}, which is marginal across the whole ${\ensuremath{M_{{\rm di}{\text -}\PJ\mskip -3mu/\mskip -2mu\Ppsi}}}\xspace$ range. The efficiency is estimated in the same way as described in Ref.~\cite{LHCb-PAPER-2016-057}.

\begin{figure}[!htbp]
\begin{center}
   \includegraphics[width=0.8\linewidth]{./figs/FigS5.pdf}
\vspace*{-0.1cm}
\end{center}
   \caption{Dependence of the experimental efficiency on ${\ensuremath{M_{{\rm di}{\text -}\PJ\mskip -3mu/\mskip -2mu\Ppsi}}}\xspace$. The error bars include the statistical uncertainties only.
   }\label{fig:eff}
\end{figure}

\clearpage

\section{Introduction}
The strong interaction is one of the fundamental forces of nature and it governs the dynamics of quarks and gluons. According to quantum chromodynamics (QCD), the theory describing the strong interaction, quarks are confined into hadrons, in agreement with experimental observations.
The quark model~\cite{GellMann:1964nj,Zweig:352337} classifies hadrons into conventional mesons~(${\ensuremath{\Pq}}\xspace{\ensuremath{\overline \quark}}\xspace$) and baryons~(${\ensuremath{\Pq}}\xspace\quark{\ensuremath{\Pq}}\xspace$ or ${\ensuremath{\overline \quark}}\xspace\quarkbar{\ensuremath{\overline \quark}}\xspace$), and also allows for the existence of
exotic hadrons such as tetraquarks~(${\ensuremath{\Pq}}\xspace\quark{\ensuremath{\overline \quark}}\xspace\quarkbar$) and pentaquarks~(${\ensuremath{\Pq}}\xspace\quark{\ensuremath{\Pq}}\xspace{\ensuremath{\overline \quark}}\xspace\quarkbar$). Exotic states provide a unique environment to study the strong interaction and the confinement mechanism~\cite{Brambilla:2014jmp}.
The first experimental evidence for an exotic hadron candidate was the $\chi_{c1}(3872)$ state observed in 2003 by the Belle collaboration~\cite{X3872}. Since then
a series of novel states consistent with
containing four quarks have been
discovered~\cite{Olsen:2017bmm}.
Recently, the \mbox{LHCb}\xspace collaboration observed resonances interpreted to be
pentaquark states~\cite{LHCb-PAPER-2015-029,LHCb-PAPER-2016-009,LHCb-PAPER-2016-015,LHCb-PAPER-2019-014}.
All hadrons observed to date, including those of exotic nature,
contain at most two heavy charm~(${\ensuremath{\Pc}}\xspace$) or bottom~(${\ensuremath{\Pb}}\xspace$) quarks, whereas many QCD-motivated phenomenological models also predict
the existence of states consisting of four heavy quarks,
\mbox{\itshape i.e.}\xspace $T_{Q_1 Q_2 \overline{Q}_3 \overline{Q}_4}$, where $Q_i$ is a {\ensuremath{\Pc}}\xspace or a {\ensuremath{\Pb}}\xspace quark~\cite{Iwasaki:1976cn,Chao:1980dv,PhysRevD.25.2370,PhysRevD.29.426,Badalian:1985es,
Berezhnoy:2011xn,Wu:2016vtq,Karliner:2016zzc,
Barnea:2006sd,Debastiani:2017msn,Liu:2019zuc,
Chen:2016jxd,
Wang:2019rdo,Bedolla:2019zwg,Lloyd:2003yc,
Chen:2020lgj,Wang:2018poa,
Anwar:2017toa,Liu:2019zoy,Esposito:2018cwh,
Becchi:2020mjz,Bai:2016int,Richard:2017vry,Vega-Morales:2017pmm,Chen:2019vrj,Berezhnoy:2012tu}.
Theoretically, the interpretation of the internal structure of such states usually assumes the formation of a $Q_1 Q_2$ diquark and a $\overline{Q}_3\overline{Q}_4$ antidiquark attracting each other.
Application of this diquark model successfully predicts
the mass of the doubly charmed
baryon ${\ensuremath{\Xires^{++}_{\cquark\cquark}}}\xspace$ \cite{Karliner:2014gca, LHCb-PAPER-2017-018} and helps to explain the relative rates of bottom baryon decays~\cite{LHCb-PAPER-2019-039}.

Tetraquark states comprising only bottom quarks, $\bbbb$, have been searched for by the \mbox{LHCb}\xspace and \mbox{CMS}\xspace collaborations
in the $\PUpsilon{\ensuremath{\Pmu^+\Pmu^-}}\xspace$ decay~\cite{LHCb-PAPER-2018-027,Sirunyan:2020txn}, with the $\PUpsilon$ state consisting of a ${\ensuremath{\Pb}}\xspace {\ensuremath{\overline \bquark}}\xspace$ pair. However, the
four-charm states, $\cccc$, have not yet been
studied in detail experimentally.
A $\cccc$ state could disintegrate into a pair of charmonium states such as ${\ensuremath{{\PJ\mskip -3mu/\mskip -2mu\Ppsi}}}\xspace$ mesons, with each consisting of a ${\ensuremath{\Pc}}\xspace {\ensuremath{\overline \cquark}}\xspace$ pair. Decays to
a {\ensuremath{{\PJ\mskip -3mu/\mskip -2mu\Ppsi}}}\xspace meson plus a heavier charmonium state, or two heavier charmonium states, with the heavier states decaying subsequently into a {\ensuremath{{\PJ\mskip -3mu/\mskip -2mu\Ppsi}}}\xspace meson and accompanying particles, are also possible.
Predictions for the masses of $\cccc$ states vary from $5.8$ to $7.4\ensuremath{\aunit{Ge\kern -0.1em V\!/}c^2}\xspace$
~\cite{Iwasaki:1976cn,Chao:1980dv,PhysRevD.25.2370,PhysRevD.29.426,Badalian:1985es,Berezhnoy:2011xn,Wu:2016vtq,Karliner:2016zzc,Barnea:2006sd,Debastiani:2017msn,Liu:2019zuc,Chen:2016jxd,Wang:2019rdo,Bedolla:2019zwg,Lloyd:2003yc,Chen:2020lgj,Wang:2018poa},
which are above the masses of
 known charmonia and charmonium-like exotic states and below those of bottomonium hadrons.~\footnote{Energy units $\aunit{Me\kern -0.1em V}\xspace=10^{6}\aunit{e\kern -0.1em V}\xspace$, $\aunit{Ge\kern -0.1em V}\xspace=10^{9}\aunit{e\kern -0.1em V}\xspace$ and $\aunit{Te\kern -0.1em V}\xspace=10^{12}\aunit{e\kern -0.1em V}\xspace$ are used in this paper.}
This mass range guarantees a clean experimental environment to identify possible $\cccc$ states
in the ${\ensuremath{{\PJ\mskip -3mu/\mskip -2mu\Ppsi}}}\xspace$-pair (also referred to as di-${\ensuremath{{\PJ\mskip -3mu/\mskip -2mu\Ppsi}}}\xspace$) invariant mass~(${\ensuremath{M_{{\rm di}{\text -}\PJ\mskip -3mu/\mskip -2mu\Ppsi}}}\xspace$) spectrum.

In proton-proton ($pp$) collisions, a pair of ${\ensuremath{{\PJ\mskip -3mu/\mskip -2mu\Ppsi}}}\xspace$ mesons can be produced in two separate interactions of gluons or quarks, named double-parton scattering\,(DPS)~\cite{Calucci:1997ii,Calucci:1999yz,DelFabbro:2000ds}, or in a single interaction, named single-parton scattering\,(SPS)~\cite{Sun:2014gca,Likhoded:2016zmk,Shao:2012iz,Shao:2015vga,Baranov:2011zz,Lansberg:2013qka,Lansberg:2014swa,Lansberg:2015lva}. The SPS process
includes both resonant production via intermediate states, which could be $\cccc$ tetraquarks, and nonresonant production.
Within the DPS process,
the two ${\ensuremath{{\PJ\mskip -3mu/\mskip -2mu\Ppsi}}}\xspace$ mesons are usually assumed to be produced independently,
thus the distribution of any di-${\ensuremath{{\PJ\mskip -3mu/\mskip -2mu\Ppsi}}}\xspace$ observable
can be constructed using the kinematics from single ${\ensuremath{{\PJ\mskip -3mu/\mskip -2mu\Ppsi}}}\xspace$ production.
Evidence of DPS has been found in
studies at the Large Hadron Collider (\mbox{LHC}\xspace) experiments~\cite{
W2jets,CMSDJ,DUpsilon,JpsiW,JpsiZ,Aaboud:2016dea,ATLASDJ,LHCb-PAPER-2012-003,LHCb-PAPER-2013-062,LHCb-PAPER-2015-046} and the AFS experiment~\cite{4jets1} in $pp$ collisions,
and at the Tevatron\xspace experiments~\cite{Abe:1993rv,Gamma3jets1,Gamma3jets2,Abazov:2014fha,Abazov:2015nnn} and the UA2 experiment~\cite{4jets2} in proton-antiproton collisions.
The \mbox{LHCb}\xspace experiment has measured
the di-${\ensuremath{{\PJ\mskip -3mu/\mskip -2mu\Ppsi}}}\xspace$ production in $pp$ collisions
at centre-of-mass energies of $\ensuremath{\protect\sqrt{s}}\xspace=7$~\cite{LHCb-PAPER-2011-013} and $13\aunit{Te\kern -0.1em V}\xspace$~\cite{LHCb-PAPER-2016-057}.
The DPS contribution is found to dominate the high
${\ensuremath{M_{{\rm di}{\text -}\PJ\mskip -3mu/\mskip -2mu\Ppsi}}}\xspace$ region, in agreement with expectation.

In this paper, fully charmed
tetraquark states $\cccc$ are searched for in the di-${\ensuremath{{\PJ\mskip -3mu/\mskip -2mu\Ppsi}}}\xspace$ invariant mass spectrum,
using $pp$ collision data collected by \mbox{LHCb}\xspace
at $\ensuremath{\protect\sqrt{s}}\xspace=7, 8$
and $13\aunit{Te\kern -0.1em V}\xspace$,
corresponding to an integrated luminosity of $9\ensuremath{\fb^{-1}}\xspace$.
The two ${\ensuremath{{\PJ\mskip -3mu/\mskip -2mu\Ppsi}}}\xspace$ candidates in a pair are reconstructed through the ${\ensuremath{{\PJ\mskip -3mu/\mskip -2mu\Ppsi}}}\xspace\ensuremath{\rightarrow}\xspace{\ensuremath{\Pmu^+}}\xspace\mun$ decay, and are labelled randomly as either ${\ensuremath{{\PJ\mskip -3mu/\mskip -2mu\Ppsi_1\mskip 2mu}}}\xspace$ or ${\ensuremath{{\PJ\mskip -3mu/\mskip -2mu\Ppsi_2\mskip 2mu}}}\xspace$.

\section{Detector and data set}
The LHCb detector is designed to study particles containing {\ensuremath{\Pb}}\xspace\ or {\ensuremath{\Pc}}\xspace\ quarks
at the \mbox{LHC}\xspace.
It is a single-arm forward spectrometer covering the pseudorapidity range $2 < \eta < 5$,
described in detail in Refs.~\cite{LHCb-DP-2008-001,LHCb-DP-2014-002}.
The online event selection is performed by a trigger,
which consists of a hardware stage, based on information from the calorimeter and muon
systems, followed by a software stage, which applies a full event
reconstruction.
At the hardware stage,
events are required to have at least one muon with high momentum transverse to the beamline, $\ensuremath{p_{\mathrm{T}}}\xspace$.
At the software stage,
two oppositely charged muon candidates
are required to have high \ensuremath{p_{\mathrm{T}}}\xspace and to form a common vertex.
Events are retained if there is at least one ${\ensuremath{{\PJ\mskip -3mu/\mskip -2mu\Ppsi}}}\xspace$ candidate passing both the hardware and software trigger requirements.
Imperfections in the description of the
magnetic field and misalignment of subdetectors lead to a bias in the momentum
measurement of charged particles, which is calibrated using
reconstructed ${\ensuremath{{\PJ\mskip -3mu/\mskip -2mu\Ppsi}}}\xspace$ and $B^+$ mesons~\cite{LHCb-PAPER-2011-035},
with well-known masses.

Simulated ${\ensuremath{{\PJ\mskip -3mu/\mskip -2mu\Ppsi}}}\xspace\ensuremath{\rightarrow}\xspace{\ensuremath{\Pmu^+}}\xspace\mun$ decays are used
to study the signal properties, including the invariant mass resolution and the reconstruction efficiency.
In the simulation, $pp$ collisions are generated using
\mbox{\textsc{Pythia}}\xspace~\cite{Sjostrand:2007gs}
 with a specific \mbox{LHCb}\xspace
configuration~\cite{LHCb-PROC-2010-056}.
Decays of unstable particles
are described by \mbox{\textsc{EvtGen}}\xspace~\cite{Lange:2001uf}, in which final-state
radiation is generated using \mbox{\textsc{Photos}}\xspace~\cite{Golonka:2005pn}.
The interaction of the generated particles with the detector and its response
are implemented using the \mbox{\textsc{Geant4}}\xspace
toolkit~\cite{Allison:2006ve}, as described in
Ref.~\cite{LHCb-PROC-2011-006}.

\section{Candidate selection}
In the offline selection,
two pairs of oppositely charged muon candidate tracks
are reconstructed,
with each pair forming a vertex of a ${\ensuremath{{\PJ\mskip -3mu/\mskip -2mu\Ppsi}}}\xspace$ candidate.
Each muon track must have \mbox{$\ensuremath{p_{\mathrm{T}}}\xspace>0.65\ensuremath{\aunit{Ge\kern -0.1em V\!/}c}\xspace$} and momentum \mbox{$p>6\ensuremath{\aunit{Ge\kern -0.1em V\!/}c}\xspace$}.
The ${\ensuremath{{\PJ\mskip -3mu/\mskip -2mu\Ppsi}}}\xspace$ candidates are required to have a dimuon invariant mass in the range $3.0<M_{\mu\mu}<3.2\ensuremath{\aunit{Ge\kern -0.1em V\!/}c^2}\xspace$.
A kinematic fit is performed for each {\ensuremath{{\PJ\mskip -3mu/\mskip -2mu\Ppsi}}}\xspace candidate constraining its vertex to coincide with a primary  $pp$ collision vertex (PV)~\cite{Hulsbergen:2005pu}.
The requirement of a good kinematic fit quality strongly suppresses
the contamination of di-${\ensuremath{{\PJ\mskip -3mu/\mskip -2mu\Ppsi}}}\xspace$ candidates stemming from feed-down of ${\ensuremath{\Pb}}\xspace$-hadrons,
which decay at displaced vertices.
The four muon tracks in  a ${\ensuremath{{\PJ\mskip -3mu/\mskip -2mu\Ppsi}}}\xspace$-pair candidate are required to originate from the same PV, reducing to a negligible level the number of pile-up candidates with the two ${\ensuremath{{\PJ\mskip -3mu/\mskip -2mu\Ppsi}}}\xspace$ candidates produced in
separated $pp$ collisions.
Fake di-${\ensuremath{{\PJ\mskip -3mu/\mskip -2mu\Ppsi}}}\xspace$ candidates, comprising two muon-track candidates reconstructed from the same real particle, are rejected by requiring muons of the same charge to have trajectories separated by an angle inconsistent with zero.
For events with more than one reconstructed di-${\ensuremath{{\PJ\mskip -3mu/\mskip -2mu\Ppsi}}}\xspace$ candidate, accounting for about 0.8\% of the total sample,
only one pair is randomly chosen.

\begin{figure}[t]
\begin{center}
\includegraphics[width=1.0\linewidth]{./figs/Fig1.pdf}
\vspace*{-1.cm}
\end{center}
   \caption{\small
	(Bottom right) Two-dimensional \mbox{$(M^{(1)}_{\mu\mu},M^{(2)}_{\mu\mu})$} distribution
   of di-${\ensuremath{{\PJ\mskip -3mu/\mskip -2mu\Ppsi}}}\xspace$ candidates
   and its projections on (bottom left)~$M^{(1)}_{\mu\mu}$ and (top)~$M^{(2)}_{\mu\mu}$. Four components are present as each projection consists of signal and background ${\ensuremath{{\PJ\mskip -3mu/\mskip -2mu\Ppsi}}}\xspace$ candidates. The labels ${\ensuremath{{\PJ\mskip -3mu/\mskip -2mu\Ppsi}}}\xspace_{1,2}$ and bkg${}_{1,2}$ represent the signal and background contributions, respectively, in the $M^{(1),(2)}_{\mu\mu}$ distribution.
   }\label{fig:2DM}
\end{figure}

The di-${\ensuremath{{\PJ\mskip -3mu/\mskip -2mu\Ppsi}}}\xspace$ signal yield is extracted by
performing an extended unbinned maximum-likelihood fit to
the two\nobreakdash-dimensional distribution of ${\ensuremath{{\PJ\mskip -3mu/\mskip -2mu\Ppsi_1\mskip 2mu}}}\xspace$ and ${\ensuremath{{\PJ\mskip -3mu/\mskip -2mu\Ppsi_2\mskip 2mu}}}\xspace$ invariant masses, \mbox{$(M^{(1)}_{\mu\mu},M^{(2)}_{\mu\mu})$},
as displayed in Fig.~\ref{fig:2DM},
where projections of the data and the fit result are shown.
For both ${\ensuremath{{\PJ\mskip -3mu/\mskip -2mu\Ppsi}}}\xspace$ candidates,
the signal mass shape is modelled by a Gaussian kernel with power-law tails~\cite{Skwarnicki:1986xj}.
Each component of combinatorial background,
consisting of random combinations of muon tracks, is described by an exponential function.
The total di-${\ensuremath{{\PJ\mskip -3mu/\mskip -2mu\Ppsi}}}\xspace$ signal yield is measured to be \mbox{$(33.57\pm0.23)\times10^3$}, where the uncertainty is statistical.

The di-${\ensuremath{{\PJ\mskip -3mu/\mskip -2mu\Ppsi}}}\xspace$ transverse momentum (${\ensuremath{p^{{\rm di}{\text -}\PJ\mskip -3mu/\mskip -2mu\Ppsi}_{\rm T}}}\xspace$) in SPS production is expected to be, on average, higher than that in DPS~\cite{Lansberg:2014swa}.
The high-${\ensuremath{p^{{\rm di}{\text -}\PJ\mskip -3mu/\mskip -2mu\Ppsi}_{\rm T}}}\xspace$ region is thus exploited to select a data sample with enhanced SPS production, which could include contributions from $\cccc$ states. Two different approaches are applied.
In the first approach~(denoted as ${\ensuremath{p^{{\rm di}{\text -}\PJ\mskip -3mu/\mskip -2mu\Ppsi}_{\rm T}}}\xspace$-threshold), ${\ensuremath{{\PJ\mskip -3mu/\mskip -2mu\Ppsi}}}\xspace$-pair candidates are selected with the requirement  ${\ensuremath{p^{{\rm di}{\text -}\PJ\mskip -3mu/\mskip -2mu\Ppsi}_{\rm T}}}\xspace>5.2\ensuremath{\aunit{Ge\kern -0.1em V\!/}c}\xspace$,
which maximises the statistical significance of the SPS signal yield, $N_{\rm SPS}/\sqrt{N_{\rm total}}$.
$N_{\rm SPS}$ and $N_{\rm total}$ are yields of the SPS component and total di-${\ensuremath{{\PJ\mskip -3mu/\mskip -2mu\Ppsi}}}\xspace$ candidates in the ${\ensuremath{M_{{\rm di}{\text -}\PJ\mskip -3mu/\mskip -2mu\Ppsi}}}\xspace$ range between $6.2$ and $7.4\ensuremath{\aunit{Ge\kern -0.1em V\!/}c^2}\xspace$, respectively. This mass region covers the predicted masses of  $\cccc$ states decaying into a ${\ensuremath{{\PJ\mskip -3mu/\mskip -2mu\Ppsi}}}\xspace$ pair.
In the second approach~(denoted as ${\ensuremath{p^{{\rm di}{\text -}\PJ\mskip -3mu/\mskip -2mu\Ppsi}_{\rm T}}}\xspace$-binned), di-${\ensuremath{{\PJ\mskip -3mu/\mskip -2mu\Ppsi}}}\xspace$ candidates are categorised  into six ${\ensuremath{p^{{\rm di}{\text -}\PJ\mskip -3mu/\mskip -2mu\Ppsi}_{\rm T}}}\xspace$ intervals with boundaries $\left\{ 0, 5, 6, 8, 9.5, 12, 50 \right\}\ensuremath{\aunit{Ge\kern -0.1em V\!/}c}\xspace$,
defined to obtain equally populated bins of SPS signal events in the $6.2<{\ensuremath{M_{{\rm di}{\text -}\PJ\mskip -3mu/\mskip -2mu\Ppsi}}}\xspace<7.4\ensuremath{\aunit{Ge\kern -0.1em V\!/}c^2}\xspace$ range.
For both scenarios, the DPS yield in the $\cccc$ signal region is extrapolated from the high-${\ensuremath{M_{{\rm di}{\text -}\PJ\mskip -3mu/\mskip -2mu\Ppsi}}}\xspace$ region using the wide-range distribution
constructed from available double-differential {\ensuremath{{\PJ\mskip -3mu/\mskip -2mu\Ppsi}}}\xspace cross-sections~\cite{LHCb-PAPER-2011-003,LHCb-PAPER-2013-016,LHCb-PAPER-2015-037}
as performed in~\cite{LHCb-PAPER-2016-057}.
The high-${\ensuremath{M_{{\rm di}{\text -}\PJ\mskip -3mu/\mskip -2mu\Ppsi}}}\xspace$ region is chosen such that the SPS yield is negligible compared to DPS.
The SPS yield is obtained by subtracting the DPS contribution from the total number of ${\ensuremath{{\PJ\mskip -3mu/\mskip -2mu\Ppsi}}}\xspace$-pair signals.

The ${\ensuremath{M_{{\rm di}{\text -}\PJ\mskip -3mu/\mskip -2mu\Ppsi}}}\xspace$ distribution for candidates with \mbox{${\ensuremath{p^{{\rm di}{\text -}\PJ\mskip -3mu/\mskip -2mu\Ppsi}_{\rm T}}}\xspace>5.2\ensuremath{\aunit{Ge\kern -0.1em V\!/}c}\xspace$}
and \mbox{$3.065< M_{\mu\mu}^{(1),(2)}<3.135\ensuremath{\aunit{Ge\kern -0.1em V\!/}c^2}\xspace$} is shown in
Fig.~\ref{fig:Xmass}. The di-${\ensuremath{{\PJ\mskip -3mu/\mskip -2mu\Ppsi}}}\xspace$ mass is calculated by constraining the reconstructed mass of each ${\ensuremath{{\PJ\mskip -3mu/\mskip -2mu\Ppsi}}}\xspace$ candidate to its known value~\cite{PDG2019}.
The spectrum shows a broad structure just above twice the ${\ensuremath{{\PJ\mskip -3mu/\mskip -2mu\Ppsi}}}\xspace$ mass threshold
ranging from $6.2$ to $6.8\ensuremath{\aunit{Ge\kern -0.1em V\!/}c^2}\xspace$ (dubbed threshold enhancement in the following)
and a narrower structure at about $6.9\ensuremath{\aunit{Ge\kern -0.1em V\!/}c^2}\xspace$, referred to hereafter as $X(6900)$.
There is also a hint of another structure around $7.2\ensuremath{\aunit{Ge\kern -0.1em V\!/}c^2}\xspace$,
whereas there are no evident structures at higher invariant mass.
Several cross-checks are performed to investigate the origin of these structures and to exclude that they are experimental artifacts.
The threshold enhancement and the $X(6900)$ structure become more pronounced in higher ${\ensuremath{p^{{\rm di}{\text -}\PJ\mskip -3mu/\mskip -2mu\Ppsi}_{\rm T}}}\xspace$ intervals, and they are present
in subsamples split according to different beam or detector configurations for data collection.
The structures are not caused by the experimental  efficiency, since the efficiency variation across the whole ${\ensuremath{M_{{\rm di}{\text -}\PJ\mskip -3mu/\mskip -2mu\Ppsi}}}\xspace$ range is found to be marginal.
Residual background, in which a muon track is reused  or at least one {\ensuremath{{\PJ\mskip -3mu/\mskip -2mu\Ppsi}}}\xspace candidate is produced from a ${\ensuremath{\Pb}}\xspace$-hadron decay, is observed to have no structure.
The possible contribution of ${\ensuremath{{\PJ\mskip -3mu/\mskip -2mu\Ppsi}}}\xspace$ pairs from ${\ensuremath{\PUpsilon}}\xspace$ decays
is estimated to be negligible and distributed uniformly in the ${\ensuremath{M_{{\rm di}{\text -}\PJ\mskip -3mu/\mskip -2mu\Ppsi}}}\xspace$ distribution. In Fig.~\ref{fig:Xmass}, the ${\ensuremath{M_{{\rm di}{\text -}\PJ\mskip -3mu/\mskip -2mu\Ppsi}}}\xspace$ distribution for background pairs with
$M_{\mu\mu}^{(1),(2)}$ in the range $3.00-3.05\ensuremath{\aunit{Ge\kern -0.1em V\!/}c^2}\xspace$ or $3.15-3.20\ensuremath{\aunit{Ge\kern -0.1em V\!/}c^2}\xspace$
is also shown, with the yield
normalised by interpolating the background into the ${\ensuremath{{\PJ\mskip -3mu/\mskip -2mu\Ppsi}}}\xspace$ signal region,
which accounts for around 15\% of the total candidates.
There is no evidence of structures in the ${\ensuremath{M_{{\rm di}{\text -}\PJ\mskip -3mu/\mskip -2mu\Ppsi}}}\xspace$ distribution of background candidates.

\begin{figure}[htb!]
\begin{center}
\includegraphics[width=0.6\linewidth]{./figs/Fig2.pdf}
\vspace*{-0.1cm}
\end{center}
   \caption{\small Invariant mass spectrum of ${\ensuremath{{\PJ\mskip -3mu/\mskip -2mu\Ppsi}}}\xspace$-pair  candidates passing the ${\ensuremath{p^{{\rm di}{\text -}\PJ\mskip -3mu/\mskip -2mu\Ppsi}_{\rm T}}}\xspace>5.2\ensuremath{\aunit{Ge\kern -0.1em V\!/}c}\xspace$ requirement with reconstructed ${\ensuremath{{\PJ\mskip -3mu/\mskip -2mu\Ppsi}}}\xspace$ masses in the  (black)~signal and (blue)~background regions,
   respectively.
   }\label{fig:Xmass}
\end{figure}

\section{Investigation of the ${\ensuremath{{\PJ\mskip -3mu/\mskip -2mu\Ppsi}}}\xspace$-pair invariant mass spectrum}
To remove background pairs that have at least one background {\ensuremath{{\PJ\mskip -3mu/\mskip -2mu\Ppsi}}}\xspace candidate, the \mbox{\em sPlot}\xspace weighting method~\cite{Pivk:2004ty} is applied, where the weights are calculated from the fit to the two\nobreakdash-dimensional  \mbox{$(M^{(1)}_{\mu\mu},M^{(2)}_{\mu\mu})$} distribution.
The background-subtracted di-${\ensuremath{{\PJ\mskip -3mu/\mskip -2mu\Ppsi}}}\xspace$ spectra in the range \mbox{$6.2<{\ensuremath{M_{{\rm di}{\text -}\PJ\mskip -3mu/\mskip -2mu\Ppsi}}}\xspace<9.0\ensuremath{\aunit{Ge\kern -0.1em V\!/}c^2}\xspace$} are shown in Fig.~\ref{fig:fitXcut} for candidates with ${\ensuremath{p^{{\rm di}{\text -}\PJ\mskip -3mu/\mskip -2mu\Ppsi}_{\rm T}}}\xspace>5.2\ensuremath{\aunit{Ge\kern -0.1em V\!/}c}\xspace$ and~Fig.~\ref{fig:fitXbinned_2peak} for candidates in the six ${\ensuremath{p^{{\rm di}{\text -}\PJ\mskip -3mu/\mskip -2mu\Ppsi}_{\rm T}}}\xspace$ intervals,
which are investigated by weighted unbinned maximum-likelihood fits~\cite{Xie:2009rka}.
The ${\ensuremath{M_{{\rm di}{\text -}\PJ\mskip -3mu/\mskip -2mu\Ppsi}}}\xspace$ distribution of signal events is expected to be dominated by the sum of the nonresonant SPS~(NRSPS) and DPS production, which have smooth shapes~(referred to as continuum in the following).
The DPS continuum is described by
a two-body phase-space function multiplied by the product of an exponential function and a second order polynomial function, whose parameters are fixed according to the ${\ensuremath{M_{{\rm di}{\text -}\PJ\mskip -3mu/\mskip -2mu\Ppsi}}}\xspace$ distribution constructed from ${\ensuremath{{\PJ\mskip -3mu/\mskip -2mu\Ppsi}}}\xspace$ differential cross-sections. Its yield is determined by extrapolation from the ${\ensuremath{M_{{\rm di}{\text -}\PJ\mskip -3mu/\mskip -2mu\Ppsi}}}\xspace>12\ensuremath{\aunit{Ge\kern -0.1em V\!/}c^2}\xspace$ region, which is dominated and well described by the DPS distribution.
The continuum NRSPS is modelled by
a two-body phase-space distribution multiplied by an exponential function
determined from the data.
The combination of continuum NRSPS and DPS does not provide a good description of the data, as is illustrated in Fig.~\ref{fig:fitXcut}(a).
The ${\ensuremath{M_{{\rm di}{\text -}\PJ\mskip -3mu/\mskip -2mu\Ppsi}}}\xspace$ spectrum in the data is tested against the hypothesis that only  NRSPS and DPS components are present
in the range $6.2<{\ensuremath{M_{{\rm di}{\text -}\PJ\mskip -3mu/\mskip -2mu\Ppsi}}}\xspace<7.4\ensuremath{\aunit{Ge\kern -0.1em V\!/}c^2}\xspace$ (null hypothesis) using a $\ensuremath{\chi^2}\xspace$ test statistic.
Pseudoexperiments are generated and fitted according to the null hypothesis, and the fraction of these fits with a $\ensuremath{\chi^2}\xspace$ value exceeding that in the data is converted into a significance.
Considering the sample in the ${\ensuremath{p^{{\rm di}{\text -}\PJ\mskip -3mu/\mskip -2mu\Ppsi}_{\rm T}}}\xspace>5.2\ensuremath{\aunit{Ge\kern -0.1em V\!/}c}\xspace$ region,
the null hypothesis is inconsistent with the data at $3.4$ standard deviations ($\sigma$).
A test performed simultaneously in the aforementioned six ${\ensuremath{p^{{\rm di}{\text -}\PJ\mskip -3mu/\mskip -2mu\Ppsi}_{\rm T}}}\xspace$ regions  yields a discrepancy  of $6.0\,\sigma$ with the null hypothesis.
A higher value is obtained in the latter case attributed to the presence of the structure at the same ${\ensuremath{M_{{\rm di}{\text -}\PJ\mskip -3mu/\mskip -2mu\Ppsi}}}\xspace$ location in different ${\ensuremath{p^{{\rm di}{\text -}\PJ\mskip -3mu/\mskip -2mu\Ppsi}_{\rm T}}}\xspace$ intervals.
A cross-check is performed by dividing the data into five or seven ${\ensuremath{p^{{\rm di}{\text -}\PJ\mskip -3mu/\mskip -2mu\Ppsi}_{\rm T}}}\xspace$ regions instead,
which results in significance values consistent with the nominal $6.0\,\sigma$.
The significance values are summarised in Table~\ref{tab:signif}~(Any structure beyond NRSPS plus DPS).

The structures in the ${\ensuremath{M_{{\rm di}{\text -}\PJ\mskip -3mu/\mskip -2mu\Ppsi}}}\xspace$ distribution can have various interpretations. There may be one or more resonant states $\cccc$ decaying directly into a pair of ${\ensuremath{{\PJ\mskip -3mu/\mskip -2mu\Ppsi}}}\xspace$ mesons, or $\cccc$ states decaying into a pair of ${\ensuremath{{\PJ\mskip -3mu/\mskip -2mu\Ppsi}}}\xspace$ mesons through feed-down of heavier quarkonia, for example $\cccc\ensuremath{\rightarrow}\xspace{\ensuremath{\Pchi_\cquark}}\xspace(\ensuremath{\rightarrow}\xspace{\ensuremath{{\PJ\mskip -3mu/\mskip -2mu\Ppsi}}}\xspace\gamma){\ensuremath{{\PJ\mskip -3mu/\mskip -2mu\Ppsi}}}\xspace$ where the photon escapes detection.
In the latter case, such a state would be expected to peak at a lower ${\ensuremath{M_{{\rm di}{\text -}\PJ\mskip -3mu/\mskip -2mu\Ppsi}}}\xspace$ position, close to the di-${\ensuremath{{\PJ\mskip -3mu/\mskip -2mu\Ppsi}}}\xspace$ mass threshold, and its structure would be broader compared to that from a direct decay.
This feed-down is unlikely an explanation for the narrow $X(6900)$ structure.
Rescattering of two charmonium states produced by SPS close to their mass threshold may also generate a narrow structure ~\cite{Braaten:2020iye,Liu:2015fea,Xie:2016lvs,Guo:2016bjq}.
The two thresholds close to the $X(6900)$ structure could be formed by ${\ensuremath{\Pchi_{\cquark 0}}}\xspace\chiczero$ pairs at $6829.4\ensuremath{\aunit{Me\kern -0.1em V\!/}c^2}\xspace$ and ${\ensuremath{\Pchi_{\cquark 1}}}\xspace{\ensuremath{\Pchi_{\cquark 0}}}\xspace$ pairs at $6925.4\ensuremath{\aunit{Me\kern -0.1em V\!/}c^2}\xspace$, respectively.
Whereas a resonance is often described by a relativistic Breit--Wigner (BW) function~\cite{PDG2019}, the lineshape of a structure with rescattering effects taken into account is more complex. In principle, resonant production can interfere with NRSPS of the same spin-parity quantum numbers~($J^\mathrm{PC}$), resulting in a coherent sum of the two components and thus a modification of the total ${\ensuremath{M_{{\rm di}{\text -}\PJ\mskip -3mu/\mskip -2mu\Ppsi}}}\xspace$ distribution.

Two different models of the structure lineshape providing a reasonable description of the data are investigated.
The $X(6900)$ lineshape parameters and yields are derived from fits to the ${\ensuremath{p^{{\rm di}{\text -}\PJ\mskip -3mu/\mskip -2mu\Ppsi}_{\rm T}}}\xspace$-threshold sample.
Simultaneous ${\ensuremath{p^{{\rm di}{\text -}\PJ\mskip -3mu/\mskip -2mu\Ppsi}_{\rm T}}}\xspace$-binned fits
are also performed as a cross-check and the variation of lineshape parameters is considered as a source of systematic uncertainties.
Due to its low significance, the structure around $7.2\ensuremath{\aunit{Ge\kern -0.1em V\!/}c^2}\xspace$ has been neglected.

In model I,
the $X(6900)$ structure is considered as a resonance, whereas the threshold enhancement is described through
a superposition of two resonances.
The lineshapes of these resonances are described by $S$-wave relativistic BW functions
multiplied by a two-body phase-space distribution.
The experimental resolution on ${\ensuremath{M_{{\rm di}{\text -}\PJ\mskip -3mu/\mskip -2mu\Ppsi}}}\xspace$ is below $5\ensuremath{\aunit{Me\kern -0.1em V\!/}c^2}\xspace$ over the full mass range and negligible compared to the widths of the structures.
The projections of the ${\ensuremath{p^{{\rm di}{\text -}\PJ\mskip -3mu/\mskip -2mu\Ppsi}_{\rm T}}}\xspace$-threshold fit using this model
are shown in
Fig.~\ref{fig:fitXcut}(b).
The mass, natural width and yield
are determined to be
\mbox{$m[X(6900)]=6905\pm11\ensuremath{\aunit{Me\kern -0.1em V\!/}c^2}\xspace$}, \mbox{$\Gamma[X(6900)]=80\pm19\aunit{Me\kern -0.1em V}\xspace$} and \mbox{$N_{\rm sig}=252\pm63$}, where biases on the statistical uncertainties have been corrected using a bootstrap method~\cite{Langenbruch:2019nwe}.
The goodness of fit is studied using a $\ensuremath{\chi^2}\xspace$ test statistic and found to be $\ensuremath{\chi^2}\xspace/{\rm ndof}=112.7/89$,
corresponding to a probability of $4.6\%$.
The fit is also performed assuming the threshold enhancement as due to  a single wide resonance (see Supplementary Material); the fit quality is found  significantly poorer and thus this model is not further investigated.

\begin{figure}[htb!]
\begin{center}
   \sidesubfloat[]{\includegraphics[width=0.48\linewidth]{./figs/Fig3a.pdf}}\\
   \sidesubfloat[]{\includegraphics[width=0.48\linewidth]{./figs/Fig3b.pdf}}
   \sidesubfloat[]{\includegraphics[width=0.48\linewidth]{./figs/Fig3c.pdf}}
\vspace*{-0.1cm}
\end{center}
   \caption{\small Invariant mass spectra of weighted di-${\ensuremath{{\PJ\mskip -3mu/\mskip -2mu\Ppsi}}}\xspace$ candidates with ${\ensuremath{p^{{\rm di}{\text -}\PJ\mskip -3mu/\mskip -2mu\Ppsi}_{\rm T}}}\xspace>5.2\ensuremath{\aunit{Ge\kern -0.1em V\!/}c}\xspace$ and overlaid projections of the ${\ensuremath{p^{{\rm di}{\text -}\PJ\mskip -3mu/\mskip -2mu\Ppsi}_{\rm T}}}\xspace$-threshold fit using (a)~the NRSPS plus DPS model, (b) model I,
   and (c) model II.
   }\label{fig:fitXcut}
\end{figure}

A comparison between the best fit result of model I
and the data reveals a tension around $6.75\ensuremath{\aunit{Ge\kern -0.1em V\!/}c^2}\xspace$, where the data shows a dip.
In an attempt to describe the dip,
model II
allows for interference between the NRSPS component and a resonance for the threshold enhancement. The coherent sum of the two components is defined as
\begin{equation}
   \label{eq:itf}
   \left| A e^{i\phi} \sqrt{f_{\rm nr}({\ensuremath{M_{{\rm di}{\text -}\PJ\mskip -3mu/\mskip -2mu\Ppsi}}}\xspace)} + {\rm BW}({\ensuremath{M_{{\rm di}{\text -}\PJ\mskip -3mu/\mskip -2mu\Ppsi}}}\xspace) \right|^2,
\end{equation}
where $A$ and $\phi$ are the magnitude and phase of the nonresonant component, relative to the BW lineshape  for the  resonance, assumed to be independent of ${\ensuremath{M_{{\rm di}{\text -}\PJ\mskip -3mu/\mskip -2mu\Ppsi}}}\xspace$,
and $f_{\rm nr}({\ensuremath{M_{{\rm di}{\text -}\PJ\mskip -3mu/\mskip -2mu\Ppsi}}}\xspace)$ is an exponential function.
The interference term in Eq.~(\ref{eq:itf}) is then added incoherently to the BW function describing the X(6900) structure and the DPS description.
The fit to the ${\ensuremath{p^{{\rm di}{\text -}\PJ\mskip -3mu/\mskip -2mu\Ppsi}_{\rm T}}}\xspace$-threshold sample with this model has a probability of $15.5\%$~($\ensuremath{\chi^2/\mathrm{ndf}}\xspace=104.7/91$), and
its projections are illustrated in Fig.~\ref{fig:fitXcut}(c).
In this case, the mass, natural width and yield
are determined to be
\mbox{$m[X(6900)]=6886\pm11\ensuremath{\aunit{Me\kern -0.1em V\!/}c^2}\xspace$}, \mbox{$\Gamma[X(6900)]=168\pm33\aunit{Me\kern -0.1em V}\xspace$} and \mbox{$N_{\rm sig}=784\pm148$}.
A larger $X(6900)$ width and yield are preferred in comparison to
model I.
Here it is assumed that the whole NRSPS production is involved in the interference with the lower-mass resonance despite that there may be several  components with different quantum numbers in the NRSPS and more than one resonance in the threshold enhancement.

Fits to the ${\ensuremath{M_{{\rm di}{\text -}\PJ\mskip -3mu/\mskip -2mu\Ppsi}}}\xspace$ distributions in the six individual ${\ensuremath{p^{{\rm di}{\text -}\PJ\mskip -3mu/\mskip -2mu\Ppsi}_{\rm T}}}\xspace$ bins are shown in Fig.~\ref{fig:fitXbinned_2peak} for
model I,
while those for
model II
are given in the Supplementary Material.
An additional model describing the dip and the $X(6900)$ structure simultaneously by using the interference between the NRSPS and a BW resonance around $6.9\ensuremath{\aunit{Ge\kern -0.1em V\!/}c^2}\xspace$ is also considered, however the fit quality is clearly poorer, as illustrated in the Supplementary Material.
Alternative lineshapes, other than the BW, may also be possible to describe these structures and will be the subject of future studies.

\begin{figure}[!tb]
\begin{center}
   \sidesubfloat[]{\includegraphics[width=0.45\linewidth]{./figs/Fig4a.pdf}}
   \sidesubfloat[]{\includegraphics[width=0.45\linewidth]{./figs/Fig4b.pdf}}\\
   \sidesubfloat[]{\includegraphics[width=0.45\linewidth]{./figs/Fig4c.pdf}}
   \sidesubfloat[]{\includegraphics[width=0.45\linewidth]{./figs/Fig4d.pdf}}\\
   \sidesubfloat[]{\includegraphics[width=0.45\linewidth]{./figs/Fig4e.pdf}}
   \sidesubfloat[]{\includegraphics[width=0.45\linewidth]{./figs/Fig4f.pdf}}\\
\vspace*{-0.1cm}
\end{center}
   \caption{\small Invariant mass spectra of weighted di-${\ensuremath{{\PJ\mskip -3mu/\mskip -2mu\Ppsi}}}\xspace$ candidates in bins of ${\ensuremath{p^{{\rm di}{\text -}\PJ\mskip -3mu/\mskip -2mu\Ppsi}_{\rm T}}}\xspace$ and overlaid projections of the ${\ensuremath{p^{{\rm di}{\text -}\PJ\mskip -3mu/\mskip -2mu\Ppsi}_{\rm T}}}\xspace$-binned fit with
  model I.
   }\label{fig:fitXbinned_2peak}
\end{figure}

The increase of the likelihood between the fits with or without considering the $X(6900)$ and the threshold enhancement structures on top of the continuum NRSPS plus DPS model is taken as the test statistic to calculate the combined global significance of the two structures~\cite{Cowan:2010js}
in the \mbox{$6.2<{\ensuremath{M_{{\rm di}{\text -}\PJ\mskip -3mu/\mskip -2mu\Ppsi}}}\xspace<7.4\ensuremath{\aunit{Ge\kern -0.1em V\!/}c^2}\xspace$} region,
where pseudoexperiments are also generated to evaluate the significance. Only model I is studied, where the interference between the NRSPS
and the threshold enhancement
is not included.
Similarly, the significance for either the threshold enhancement or the $X(6900)$ structure is evaluated assuming the presence of the other along with the NRSPS and DPS continuum.
The significance is determined from
both ${\ensuremath{p^{{\rm di}{\text -}\PJ\mskip -3mu/\mskip -2mu\Ppsi}_{\rm T}}}\xspace$-threshold and ${\ensuremath{p^{{\rm di}{\text -}\PJ\mskip -3mu/\mskip -2mu\Ppsi}_{\rm T}}}\xspace$-binned fits,
and summarised in Table~\ref{tab:signif}. The results are above $5\,\sigma$ for the two structures, with slightly higher significance for the ${\ensuremath{p^{{\rm di}{\text -}\PJ\mskip -3mu/\mskip -2mu\Ppsi}_{\rm T}}}\xspace$-binned case.

\begin{table}[!tb]
\caption{\small Global significance evaluated under the various assumptions described in the text.
   }\label{tab:signif}
\begin{center}
\begin{tabular}{c|c|c}
\toprule
   \multirow{2}{*}{Structure} & \multicolumn{2}{c}{Significance} \\
   & ${\ensuremath{p^{{\rm di}{\text -}\PJ\mskip -3mu/\mskip -2mu\Ppsi}_{\rm T}}}\xspace$-threshold & ${\ensuremath{p^{{\rm di}{\text -}\PJ\mskip -3mu/\mskip -2mu\Ppsi}_{\rm T}}}\xspace$-binned \\
\hline
   Any structure beyond NRSPS plus DPS & $3.4\,\sigma$ & $6.0\,\sigma$ \\
\hline
   Threshold enhancement plus $X(6900)$ & $6.4\,\sigma$ & $6.9\,\sigma$ \\
   Threshold enhancement & $6.0\,\sigma$ & $6.5\,\sigma$ \\
   $X(6900)$ & $5.1\,\sigma$ & $5.4\,\sigma$ \\
\bottomrule
\end{tabular}
\end{center}
\end{table}

\begin{table}[!tb]
   \caption{\small Systematic uncertainties on the mass~($m$) and natural width~($\Gamma$) of the $X(6900)$ structure.
   }\label{tab:sys}
\begin{center}
\begin{tabular}{c|cc|cc}
\toprule
   & \multicolumn{2}{c|}{Without interference} & \multicolumn{2}{c}{With interference} \\
   Component & $m$~[$\ensuremath{\aunit{Me\kern -0.1em V\!/}c^2}\xspace$] & $\Gamma$~[$\aunit{Me\kern -0.1em V}\xspace$] & $m$~[$\ensuremath{\aunit{Me\kern -0.1em V\!/}c^2}\xspace$] & $\Gamma$~[$\aunit{Me\kern -0.1em V}\xspace$] \\
\hline
\mbox{\em sPlot}\xspace weights             & 0.8 & 10.3 & 4.4 & 36.9 \\
Experimental resolution   & 0.0 & 1.4  & 0.0 & 0.6 \\
NRSPS+DPS modelling       & 0.8 & 16.1 & 3.5 & 9.3 \\
$X(6900)$ shape              & 0.0 & 0.3  & 0.4 & 0.2 \\
Dependence on ${\ensuremath{p^{{\rm di}{\text -}\PJ\mskip -3mu/\mskip -2mu\Ppsi}_{\rm T}}}\xspace$       & 4.6 & 13.5 & 6.2 & 56.7 \\
${\ensuremath{\Pb}}\xspace$-hadron feed-down & 0.0 & 0.2  & 0.0 & 5.3 \\
Structure at $7.2\ensuremath{\aunit{Ge\kern -0.1em V\!/}c^2}\xspace$          & 1.3 & 9.2  & 6.7 & 5.2 \\
Threshold structure shape & 5.2 & 20.5 & --  & -- \\
NRSPS phase               & --  & --   & 0.3 & 1.3 \\
\hline
Total & 7 & 33 & 11 & 69 \\
\bottomrule
\end{tabular}
\end{center}
\end{table}

Systematic uncertainties on the measurements of the mass and natural width of the $X(6900)$ structure
are reported in Table~\ref{tab:sys}.
They include variations of the results obtained by:
including an explicit component in the ${\ensuremath{M_{{\rm di}{\text -}\PJ\mskip -3mu/\mskip -2mu\Ppsi}}}\xspace$ fits  for the {\ensuremath{{\PJ\mskip -3mu/\mskip -2mu\Ppsi}}}\xspace combinatorial background rather than subtracting it using the weighting method (\mbox{\em sPlot}\xspace weights in Table~\ref{tab:sys});
convolving the ${\ensuremath{M_{{\rm di}{\text -}\PJ\mskip -3mu/\mskip -2mu\Ppsi}}}\xspace$ fit functions with a Gaussian function of $5\ensuremath{\aunit{Me\kern -0.1em V\!/}c^2}\xspace$ width to account for the invariant mass resolution (Experimental resolution);
using alternative functions to describe the NRSPS component  and varying the DPS yield (NRSPS plus DPS modelling);
using  an alternative $P$-wave BW function for the $X(6900)$ structure and varying the hadron radius in the BW function from $2$ to $5\aunit{Ge\kern -0.1em V}\xspace^{-1}$ [$X(6900)$ shape];
obtaining results from a simultaneous fit to the ${\ensuremath{M_{{\rm di}{\text -}\PJ\mskip -3mu/\mskip -2mu\Ppsi}}}\xspace$ distributions in the six ${\ensuremath{p^{{\rm di}{\text -}\PJ\mskip -3mu/\mskip -2mu\Ppsi}_{\rm T}}}\xspace$ bins which covers the uncertainty due to variations of the NRSPS, DPS shapes and the NRSPS-resonance interference with respect to ${\ensuremath{p^{{\rm di}{\text -}\PJ\mskip -3mu/\mskip -2mu\Ppsi}_{\rm T}}}\xspace$ (Dependence on ${\ensuremath{p^{{\rm di}{\text -}\PJ\mskip -3mu/\mskip -2mu\Ppsi}_{\rm T}}}\xspace$);
including an explicit contribution for {\ensuremath{{\PJ\mskip -3mu/\mskip -2mu\Ppsi}}}\xspace mesons from ${\ensuremath{\Pb}}\xspace$-hadron feed-down ($b$-hadron feed-down) or adding a BW component for the $7.2\ensuremath{\aunit{Ge\kern -0.1em V\!/}c^2}\xspace$ structure (Structure at $7.2\ensuremath{\aunit{Ge\kern -0.1em V\!/}c^2}\xspace$);
modelling the threshold structure using an alternative Gaussian function with asymmetric power-law tails,
or fitting in a reduced ${\ensuremath{M_{{\rm di}{\text -}\PJ\mskip -3mu/\mskip -2mu\Ppsi}}}\xspace$ range excluding the threshold structure (Threshold structure shape);
allowing the relative phase in the NRSPS component to vary linearly with ${\ensuremath{M_{{\rm di}{\text -}\PJ\mskip -3mu/\mskip -2mu\Ppsi}}}\xspace$ (NRSPS phase).
The total uncertainties
are determined to be
$7\ensuremath{\aunit{Me\kern -0.1em V\!/}c^2}\xspace$ and $33\aunit{Me\kern -0.1em V}\xspace$ for the mass and natural width, respectively, without considering any interference, and
$11\ensuremath{\aunit{Me\kern -0.1em V\!/}c^2}\xspace$ and $69\aunit{Me\kern -0.1em V}\xspace$ when the interference between NRSPS and the threshold structure is introduced.

For the scenario without interference,
the production cross-section of the $X(6900)$ structure relative to that of all ${\ensuremath{{\PJ\mskip -3mu/\mskip -2mu\Ppsi}}}\xspace$ pairs~(inclusive), times the branching fraction  ${\ensuremath{\mathcal{B}}}\xspace(X(6900)\ensuremath{\rightarrow}\xspace{\ensuremath{{\PJ\mskip -3mu/\mskip -2mu\Ppsi}}}\xspace\jpsi)$, $\mathcal{R}$,
is determined in the $pp$ collision data at $\ensuremath{\protect\sqrt{s}}\xspace=13\aunit{Te\kern -0.1em V}\xspace$.
The measurement is obtained for both ${\ensuremath{{\PJ\mskip -3mu/\mskip -2mu\Ppsi}}}\xspace$ mesons in the fiducial region of
transverse momentum
below $10\ensuremath{\aunit{Ge\kern -0.1em V\!/}c}\xspace$
and rapidity
between $2.0$ and $4.5$.
An event-by-event efficiency correction is performed to obtain the signal yield at production.
The residual contamination from ${\ensuremath{\Pb}}\xspace$-hadron feed-down is subtracted from inclusive ${\ensuremath{{\PJ\mskip -3mu/\mskip -2mu\Ppsi}}}\xspace$-pair production following
Ref.~\cite{LHCb-PAPER-2015-037}.
The systematic uncertainties on the $X(6900)$ yield
are estimated in a similar way to that for the mass and natural width,
while other systematic uncertainties mostly cancel in the ratio.
The production ratio is measured to be
\mbox{$\mathcal{R}=[1.1\pm0.4\aunit{(stat)}\xspace\pm0.3\aunit{(syst)}\xspace]\%$} without any ${\ensuremath{p^{{\rm di}{\text -}\PJ\mskip -3mu/\mskip -2mu\Ppsi}_{\rm T}}}\xspace$ requirement and \mbox{$\mathcal{R}=[2.6\pm0.6\aunit{(stat)}\xspace\pm0.8\aunit{(syst)}\xspace]\%$} for ${\ensuremath{p^{{\rm di}{\text -}\PJ\mskip -3mu/\mskip -2mu\Ppsi}_{\rm T}}}\xspace>5.2\ensuremath{\aunit{Ge\kern -0.1em V\!/}c}\xspace$.

\section{Summary}
In conclusion,
using $pp$ collision data
at centre-of-mass energies of $7$, $8$ and $13\aunit{Te\kern -0.1em V}\xspace$ collected with the \mbox{LHCb}\xspace detector,
corresponding to an integrated luminosity of $9\ensuremath{\fb^{-1}}\xspace$,
the ${\ensuremath{{\PJ\mskip -3mu/\mskip -2mu\Ppsi}}}\xspace$-pair invariant mass spectrum is studied.
The data in the mass range between $6.2$ and $7.4\ensuremath{\aunit{Ge\kern -0.1em V\!/}c^2}\xspace$ are found to be inconsistent with the hypothesis of
NRSPS plus DPS continuum.
A narrow structure, $X(6900)$, matching the lineshape of a resonance
and a broad structure next to the di-{\ensuremath{{\PJ\mskip -3mu/\mskip -2mu\Ppsi}}}\xspace mass threshold
are found.
The global significance of either the broad or the $X(6900)$ structure is determined to be larger than five standard deviations.
Describing the $X(6900)$ structure with a Breit--Wigner lineshape,
its mass and natural width
are determined to be
\begin{equation*}
m[X(6900)]=6905\pm11\pm7\ensuremath{\aunit{Me\kern -0.1em V\!/}c^2}\xspace
\end{equation*}
and
\begin{equation*}
\Gamma[X(6900)]=80\pm19\pm33\aunit{Me\kern -0.1em V}\xspace,
\end{equation*}
assuming no interference with the NRSPS continuum is present,
where the first uncertainty is statistical and the second systematic.
When assuming the NRSPS continuum interferes with the broad structure close to the di-{\ensuremath{{\PJ\mskip -3mu/\mskip -2mu\Ppsi}}}\xspace mass threshold, they become
\begin{equation*}
m[X(6900)]=6886\pm11\pm11\ensuremath{\aunit{Me\kern -0.1em V\!/}c^2}\xspace
\end{equation*}
and
\begin{equation*}
\Gamma[X(6900)]=168\pm33\pm69\aunit{Me\kern -0.1em V}\xspace.
\end{equation*}
The $X(6900)$ structure could originate from a hadron state consisting of four charm quarks, $\cccc$, predicted in various tetraquark models. The broad structure close to the di-{\ensuremath{{\PJ\mskip -3mu/\mskip -2mu\Ppsi}}}\xspace mass threshold could be due to a mixture of multiple four-charm quark states or have contributions from feed-down decays of four-charm states through heavier quarkonia.
Other interpretations cannot presently be ruled out, for example the rescattering of two charmonium states produced close to their mass threshold. More data along with additional measurements, including determination of the spin-parity quantum numbers and $\ensuremath{p_{\mathrm{T}}}\xspace$ dependence of the production cross-section,
are needed to provide further information about the nature of the observed structure.