\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.







