EUROPEAN ORGANIZATION FOR NUCLEAR RESEARCH (CERN)

 CERN-EP-2021-236 LHCb-PAPER-2021-039 March 28, 2022

Searches for rare ${B}^{0}_{s}$ and ${B}^{0}$ decays into four muons

LHCb collaboration†††Authors are listed at the end of this paper.

This paper is dedicated to the memory of our friend and colleague Alexei
Vorobyev.

Searches for rare ${B}^{0}_{s}$ and ${B}^{0}$ decays into four muons are
performed using proton-proton collision data recorded by the LHCb experiment,
corresponding to an integrated luminosity of 9$\text{\,fb}^{-1}$. Direct
decays and decays via light scalar and ${J\mskip-3.0mu/\mskip-2.0mu\psi}$
resonances are considered. No evidence for the six decays searched for is
found and upper limits at the 95% confidence level on their branching
fractions ranging between $1.8\times 10^{-10}$ and $2.6\times 10^{-9}$ are
set.

Published in JHEP 03 (2022) 109

© 2024 CERN for the benefit of the LHCb collaboration. CC BY 4.0 licence.

##  1 Introduction

Decays of neutral ${B}^{0}_{s}$ and ${B}^{0}$ mesons into four muons that are
not mediated by intermediate resonances proceed through $b\rightarrow s$ and
$b\rightarrow d$ quark flavour-changing neutral current (FCNC) transitions.
These decays proceed by electroweak loop amplitudes in the Standard Model
(SM), as illustrated in Fig. 1, due to the absence of tree-level FCNCs and
hence are highly suppressed, with predicted branching fractions of
$\mathcal{B}\left({{B}^{0}_{s}}\\!\rightarrow{\mu^{+}\mu^{-}\mu^{+}\mu^{-}}\right)=\left(0.9-1.0\right)\times
10^{-10}$ and
$\mathcal{B}\left({{B}^{0}}\\!\rightarrow{\mu^{+}\mu^{-}\mu^{+}\mu^{-}}\right)=\left(0.4-4.0\right)\times
10^{-12}$ [1].

However, new particles beyond the Standard Model (BSM) may significantly
enhance these branching fractions. For example, decays via scalar and
pseudoscalar sgoldstino particles into a pair of dimuons in the Minimal
Supersymmetric Standard Model (MSSM) may lead to significant enhancements of
the branching fractions [2]. Furthermore, rare ${B}^{0}_{s}$ and ${B}^{0}$
decays into a pair of dimuons mediated by BSM light narrow scalar resonances
($a$) of the form
${{B}_{({s})}^{0}}\\!\rightarrow{a}\left({\mu^{+}\mu^{-}}\right){a}\left({\mu^{+}\mu^{-}}\right)$
naturally occur in extensions of the SM involving a new strongly interacting
sector. In particular, such models [3, 4] can account for the long-standing
tension between the SM prediction [5] and the observed [6, 7] value of the
anomalous magnetic dipole moment of the muon, as well as the widely discussed
anomalies in ${b}\\!\rightarrow{s}\ell^{+}\ell^{-}$ transitions [8, 9, 10, 11,
12, 13, 14]. This motivates a search for $B$ decays into two light scalars
with masses around 1$\text{\,Ge\kern-1.00006ptV\\!/}c^{2}$ [15].

Figure 1: Feynman diagrams for nonresonant
${{B}^{0}_{s}},{{B}^{0}}\\!\rightarrow{\mu^{+}\mu^{-}\mu^{+}\mu^{-}}$ decays.

This article presents searches for
${{B}_{({s})}^{0}}\\!\rightarrow{\mu^{+}\mu^{-}\mu^{+}\mu^{-}}$ decays using
proton-proton ($pp$) collision data recorded by the LHCb experiment in 2011,
2012 (Run 1) and 2015-2018 (Run 2) at centre-of-mass energies of 7, 8 and 13
TeV, respectively, corresponding to a total integrated luminosity of
9$\text{\,fb}^{-1}$. Previously, the most sensitive searches for
${{B}_{({s})}^{0}}\\!\rightarrow{\mu^{+}\mu^{-}\mu^{+}\mu^{-}}$ decays were
performed by the LHCb collaboration using data recorded during Run 1
corresponding to an integrated luminosity of 3$\text{\,fb}^{-1}$. No signals
were observed and the limits
$\mathcal{B}\left({{B}^{0}_{s}}\\!\rightarrow{\mu^{+}\mu^{-}\mu^{+}\mu^{-}}\right)<2.5\times
10^{-9}$ and
$\mathcal{B}\left({{B}^{0}}\\!\rightarrow{\mu^{+}\mu^{-}\mu^{+}\mu^{-}}\right)<6.9\times
10^{-10}$ were set [16]. These searches were insensitive to
${{B}_{({s})}^{0}}\\!\rightarrow{a}\left({\mu^{+}\mu^{-}}\right){a}\left({\mu^{+}\mu^{-}}\right)$
decays with $m_{a}\approx 1\text{\,Ge\kern-1.00006ptV\\!/}c^{2}$ due to
requirements imposed to remove decays via $\phi$ mesons. Now, a specific
selection is used that avoids vetoing the dimuon mass region around
$1\text{\,Ge\kern-1.00006ptV\\!/}c^{2}$.

In addition, searches are performed for
${{B}_{({s})}^{0}}\\!\rightarrow{{J\mskip-3.0mu/\mskip-2.0mu\psi}}\left({\mu^{+}\mu^{-}}\right){\mu^{+}\mu^{-}}$
decays, which are tree-level $b\rightarrow c$ transitions that proceed via $W$
boson exchange, as demonstrated in Fig. 2.

Figure 2: Feynman diagram for
${{B}^{0}_{s}},{{B}^{0}}\\!\rightarrow{{J\mskip-3.0mu/\mskip-2.0mu\psi}}\left({\mu^{+}\mu^{-}}\right){\mu^{+}\mu^{-}}$
decays.

No SM predictions for these branching fractions are currently available.
However, an estimate for the ${B}^{0}_{s}$ decay may be obtained using a
prediction for the branching fraction of the
${{B}^{0}_{s}}\rightarrow{{J\mskip-3.0mu/\mskip-2.0mu\psi}}{e^{+}e^{-}}$ decay
[17], the measured
${{J\mskip-3.0mu/\mskip-2.0mu\psi}}\rightarrow{\mu^{+}}{\mu^{-}}$ branching
fraction [18], and assuming that the
${{B}^{0}_{s}}\rightarrow{{J\mskip-3.0mu/\mskip-2.0mu\psi}}{\mu^{+}}{\mu^{-}}$
and ${{B}^{0}_{s}}\rightarrow{{J\mskip-3.0mu/\mskip-2.0mu\psi}}{e^{+}e^{-}}$
branching fractions are equal, yielding
$\mathcal{B}\left({{B}^{0}_{s}}\\!\rightarrow{{J\mskip-3.0mu/\mskip-2.0mu\psi}}\left({\mu^{+}\mu^{-}}\right){\mu^{+}\mu^{-}}\right)\sim
10^{-11}$. The equivalent ${B}^{0}$ branching fraction may be estimated by
multiplying this value by $|V_{cb}V_{cs}|/|V_{cb}V_{cd}|\approx 0.05$, leading
to
$\mathcal{B}\left({{B}^{0}}\\!\rightarrow{{J\mskip-3.0mu/\mskip-2.0mu\psi}}\left({\mu^{+}\mu^{-}}\right){\mu^{+}\mu^{-}}\right)\sim
10^{-13}$. Since these branching fractions are significantly below the
sensitivity of this analysis, the presence of a signal would indicate BSM
physics.

The branching fractions are measured using
${{B}^{0}_{s}}\\!\rightarrow{{J\mskip-3.0mu/\mskip-2.0mu\psi}}\left({\mu^{+}\mu^{-}}\right)\phi\left({\mu^{+}\mu^{-}}\right)$
decays as a normalisation channel, which has a measured branching fraction of
$\left(1.74\pm 0.14\right)\times 10^{-8}$ [19, 18]. Since this decay has the
same final state as the signal modes, many sources of systematic uncertainty
cancel. First, the yield of the
${{B}^{0}_{s}}\\!\rightarrow{{J\mskip-3.0mu/\mskip-2.0mu\psi}}\left({\mu^{+}\mu^{-}}\right)\phi\left({\mu^{+}\mu^{-}}\right)$
control mode is estimated using an unbinned maximum-likelihood fit to its
invariant mass spectrum. Then, searches for
${{B}_{({s})}^{0}}\\!\rightarrow{\mu^{+}\mu^{-}\mu^{+}\mu^{-}}$ and
${{B}_{({s})}^{0}}\\!\rightarrow{{J\mskip-3.0mu/\mskip-2.0mu\psi}}\left({\mu^{+}\mu^{-}}\right){\mu^{+}\mu^{-}}$
decays are made using unbinned maximum-likelihood fits to the invariant mass
spectra of each decay. These are performed simultaneously in multiple
intervals of the response of a multivariate classifier constructed to separate
signal and combinatorial background, where combinatorial background arises
from muons that do not all originate from the same $b$-hadron decay. In
contrast,
${{B}_{({s})}^{0}}\\!\rightarrow{a}\left({\mu^{+}\mu^{-}}\right){a}\left({\mu^{+}\mu^{-}}\right)$
decays are searched for by imposing a minimum requirement on the multivariate
classifier response due to the very low levels of background, followed by a
single fit to the mass spectrum of each signal mode. For all six signal modes,
the branching fractions are normalised using the yield of the control mode,
its measured branching fraction, and the ratio of efficiencies between the
signal and control modes calculated using simulated events. In the case of
${B}^{0}$ decays, the ratio of fragmentation fractions $f_{s}/f_{d}$ is also
taken into account.

##  2 Detector and simulation

The LHCb detector is a single-arm forward spectrometer covering the
pseudorapidity range $2<\eta<5$, described in detail in Refs. [20, 21]. It
includes a high-precision tracking system consisting of a silicon-strip vertex
locator (VELO), surrounding the $pp$ interaction region, a large-area silicon-
strip detector located upstream of a dipole magnet with a bending power of
about 4 Tm, and three stations of silicon-strip detectors and straw drift
tubes placed downstream of the magnet. Particle identification is provided by
two ring-imaging Cherenkov (RICH) detectors, an electromagnetic and a hadronic
calorimeter, and a muon system composed of alternating layers of iron and
multiwire proportional chambers.

In the simulation, $pp$ collisions are generated using Pythia [22,
*Sjostrand:2006za] with a specific LHCb configuration [24]. Decays of unstable
particles are described by EvtGen [25], in which final-state radiation is
generated using Photos [26]. The interaction of the generated particles with
the detector, and its response, are implemented using the Geant4 toolkit [27,
*Agostinelli:2002hh] as described in Ref. [29]. Simulated
${{B}_{({s})}^{0}}\\!\rightarrow{\mu^{+}\mu^{-}\mu^{+}\mu^{-}}$ decays are
generated using a phase-space model for the kinematics of the muons, due to
the lack of a SM prediction for the decay dynamics. The
${{B}_{({s})}^{0}}\\!\rightarrow{a}\left({\mu^{+}\mu^{-}}\right){a}\left({\mu^{+}\mu^{-}}\right)$
decays are simulated with an $a$ state with
$m_{a}=1\text{\,Ge\kern-1.00006ptV\\!/}c^{2}$, a lifetime of $1\text{\,fs}$
and a natural width of zero. With these settings it is assumed that the $a$
natural width is significantly below the experimental resolution of
approximately 20$\text{\,Me\kern-1.00006ptV\\!/}c^{2}$. Simulated
${{B}_{({s})}^{0}}\\!\rightarrow{{J\mskip-3.0mu/\mskip-2.0mu\psi}}\left({\mu^{+}\mu^{-}}\right){\mu^{+}\mu^{-}}$
decays are generated using the BTOSLLBALL model [30].

##  3 Event selection

The online selection is composed of a hardware trigger, followed by a two-
stage software trigger. Candidate
${{B}_{({s})}^{0}}\\!\rightarrow{\mu^{+}\mu^{-}\mu^{+}\mu^{-}}$ decays are
required to pass triggers designed to select decays involving muons, which are
identified from tracks that penetrate the calorimeters and the iron layers
between the muon stations. In the hardware stage, candidates must pass at
least one of two selections: one requiring the presence of at least one high
transverse momentum ($p_{\mathrm{T}}$) muon, the other requiring a pair of
muons with a product of their respective transverse momenta above a threshold.
In the first software trigger stage, candidates must include a
high-$p_{\mathrm{T}}$ muon that is inconsistent with having originated at a
primary $pp$ interaction vertex (PV). In the second stage, requirements are
imposed on pairs of muons in order to select candidates consistent with $B$
meson decay vertices that are significantly displaced from the PV.

Candidate ${{B}_{({s})}^{0}}\\!\rightarrow{\mu^{+}\mu^{-}\mu^{+}\mu^{-}}$
decays are formed by combining four tracks identified as muons that originate
from a common decay vertex. The muons forming the candidate are all required
to be inconsistent with originating from a PV, have a $p_{\mathrm{T}}$ >
250$\text{\,Me\kern-1.00006ptV\\!/}c$, and a good track fit quality. The
maximum distance of closest approach between the four tracks is required to be
below 0.3 mm. The resulting $B$ candidates are required to have an invariant
mass within 1$\text{\,Ge\kern-1.00006ptV\\!/}c^{2}$ of the known ${B}^{0}_{s}$
mass [18] in the Run 1 dataset or larger than
4$\text{\,Ge\kern-1.00006ptV\\!/}c^{2}$ in the Run 2 dataset, a good quality
vertex fit, to be consistent with originating from a PV, have a significant
flight distance, and the cosine of the angle between its momentum vector and
the vector pointing from the PV to its decay vertex greater than zero.
Hadronic background is suppressed using particle identification (PID)
information provided by the muon system, calorimeters and RICH detectors,
which is used to select well identified muons [31].

The different signal modes and the control channel are separated using
requirements on the four invariant masses of pairs of oppositely charged muons
$q_{ij}$, where $i$ and $j$ index the positively and negatively charged muons,
respectively, such that $ij$ can take the values $\\{11,12,21,22\\}$. Regions
around the $\phi$, ${J\mskip-3.0mu/\mskip-2.0mu\psi}$ and $\psi{(2S)}$ mesons
are defined as
$950\text{\,Me\kern-1.00006ptV\\!/}c^{2}<q_{ij}<1090\text{\,Me\kern-1.00006ptV\\!/}c^{2}$,
$\lvert
q_{ij}-m({{J\mskip-3.0mu/\mskip-2.0mu\psi}})\rvert<100\text{\,Me\kern-1.00006ptV\\!/}c^{2}$
and $\lvert
q_{ij}-m({\psi{(2S)}})\rvert<100\text{\,Me\kern-1.00006ptV\\!/}c^{2}$,
respectively, where $m({{J\mskip-3.0mu/\mskip-2.0mu\psi}})$ and
$m({\psi{(2S)}})$ are the known ${J\mskip-3.0mu/\mskip-2.0mu\psi}$ and
$\psi{(2S)}$ masses [18]. The
${{B}_{({s})}^{0}}\\!\rightarrow{\mu^{+}\mu^{-}\mu^{+}\mu^{-}}$ signal
selection requires that none of the four pairs of opposite-sign muons has an
invariant mass within any of the $\phi$, ${J\mskip-3.0mu/\mskip-2.0mu\psi}$ or
$\psi{(2S)}$ regions. This requirement is approximately 64% efficient on
simulated ${{B}^{0}_{s}}\\!\rightarrow{\mu^{+}\mu^{-}\mu^{+}\mu^{-}}$
candidates passing the trigger and offline selection, while rejecting 99.94%
of the
${{B}^{0}_{s}}\\!\rightarrow{{J\mskip-3.0mu/\mskip-2.0mu\psi}}\left({\mu^{+}\mu^{-}}\right)\phi\left({\mu^{+}\mu^{-}}\right)$
decays. The control mode selection requires that at least one of the mutually
exclusive pairs of dimuons (i.e. $\mu_{1}^{+}\mu_{1}^{-}$ and
$\mu_{2}^{+}\mu_{2}^{-}$ or $\mu_{1}^{+}\mu_{2}^{-}$ and
$\mu_{2}^{+}\mu_{1}^{-}$) has one dimuon mass in the $\phi$ region while the
other falls in the ${J\mskip-3.0mu/\mskip-2.0mu\psi}$ region. This selection
is around 95% efficient on simulated
${{B}^{0}_{s}}\\!\rightarrow{{J\mskip-3.0mu/\mskip-2.0mu\psi}}\left({\mu^{+}\mu^{-}}\right)\phi\left({\mu^{+}\mu^{-}}\right)$
decays.

Candidate
${{B}_{({s})}^{0}}\\!\rightarrow{a}\left({\mu^{+}\mu^{-}}\right){a}\left({\mu^{+}\mu^{-}}\right)$
decays are selected by requiring two mutually exclusive pairs of opposite-sign
muons with similar invariant masses, satisfying

$\lvert
q^{2}_{ij}-q^{2}_{kl}\rvert<2\sqrt{\sigma^{2}(q^{2}_{ij})+\sigma^{2}(q^{2}_{kl})}\,,$
(1)

where $\sigma(q^{2}_{ij})$ is the dimuon invariant mass squared resolution
evaluated at $q^{2}_{ij}$. The resolution increases roughly linearly as a
function of $q^{2}$, from around 0.003 $\text{\,Ge\kern-1.00006ptV}^{2}$ to
around 0.15 $\text{\,Ge\kern-1.00006ptV}^{2}$, and is estimated as a function
of $q^{2}$ separately in each data-taking year using simulated
${{B}^{0}_{s}}\\!\rightarrow{\mu^{+}\mu^{-}\mu^{+}\mu^{-}}$ decays. This
requirement does not completely eliminate contamination from
${{B}^{0}_{s}}\\!\rightarrow{{J\mskip-3.0mu/\mskip-2.0mu\psi}}\left({\mu^{+}\mu^{-}}\right)\phi\left({\mu^{+}\mu^{-}}\right)$
decays, which may satisfy the requirement if muons from the intermediate
${J\mskip-3.0mu/\mskip-2.0mu\psi}$ and $\phi$ mesons are paired incorrectly,
giving two intermediate combinations with similar masses. This remaining
background is vetoed by requiring that none of the four pairs of oppositely
charged muons has an invariant mass in the ${J\mskip-3.0mu/\mskip-2.0mu\psi}$
region.

The
${{B}_{({s})}^{0}}\\!\rightarrow{{J\mskip-3.0mu/\mskip-2.0mu\psi}}\left({\mu^{+}\mu^{-}}\right){\mu^{+}\mu^{-}}$
candidates are selected by requiring that at least one of the pairs of
opposite-sign muons has a mass in the ${J\mskip-3.0mu/\mskip-2.0mu\psi}$
region. Meanwhile the corresponding opposite-sign pair of muons is required to
have a mass that falls outside the $\phi$ region under both a dimuon and
dikaon mass hypothesis, in order to remove background from
${{B}^{0}_{s}}\\!\rightarrow{{J\mskip-3.0mu/\mskip-2.0mu\psi}}\left({\mu^{+}\mu^{-}}\right)\phi\left({\mu^{+}\mu^{-}}\right)$
and
${{B}^{0}_{s}}\\!\rightarrow{{J\mskip-3.0mu/\mskip-2.0mu\psi}}\left({\mu^{+}\mu^{-}}\right)\phi\left({{K}^{+}}{{K}^{-}}\right)$
decays, respectively. To suppress additional hadronic background such as
${{B}^{0}_{s}}\\!\rightarrow{{J\mskip-3.0mu/\mskip-2.0mu\psi}}\left({\mu^{+}\mu^{-}}\right){{\pi}^{+}}{{\pi}^{-}}$
and
${{B}^{0}}\\!\rightarrow{{J\mskip-3.0mu/\mskip-2.0mu\psi}}\left({\mu^{+}\mu^{-}}\right){{K}^{*0}}\left({{K}^{+}}{{\pi}^{-}}\right)$
decays, the two muons not forming the ${J\mskip-3.0mu/\mskip-2.0mu\psi}$
candidate are required to satisfy more stringent PID criteria.

Finally, a multivariate classifier based on a boosted decision tree (BDT)
algorithm from the TMVA package [32] is used to separate signal from
combinatorial background. The classifier is trained using a mixture of
simulated ${{B}^{0}_{s}}\\!\rightarrow{\mu^{+}\mu^{-}\mu^{+}\mu^{-}}$ and
${{B}^{0}}\\!\rightarrow{\mu^{+}\mu^{-}\mu^{+}\mu^{-}}$ decays as the signal
proxy, while the background proxy comprises data with an invariant $B$ mass,
$m({\mu^{+}\mu^{-}}{\mu^{+}\mu^{-}})$, in the regions
$m({\mu^{+}\mu^{-}}{\mu^{+}\mu^{-}})>5426\text{\,Me\kern-1.00006ptV\\!/}c^{2}$
and
$m({\mu^{+}\mu^{-}}{\mu^{+}\mu^{-}})<5020\text{\,Me\kern-1.00006ptV\\!/}c^{2}$.
The simulated ${{B}_{({s})}^{0}}\\!\rightarrow{\mu^{+}\mu^{-}\mu^{+}\mu^{-}}$
decays used in the training are corrected to improve agreement with data as
described in Sec. 5. Two separate BDT classifiers are trained for the Run 1
and Run 2 datasets due to the different data-taking conditions. The size of
the available training sample is maximised in each case using the k-folding
technique [33], which ensures that the BDT used to classify a given event in
the data background sample is trained independent of the event itself. The
following properties of the $B$ candidate are taken as inputs to the
classifiers: the logarithm of its impact parameter $\chi^{2}$ with respect to
the associated PV, the logarithm of its flight distance $\chi^{2}$, its flight
distance, pseudorapidity, $p_{\mathrm{T}}$, decay time, decay vertex
$\chi^{2}$ per degree of freedom, and the minimum impact parameter $\chi^{2}$
with respect to the associated PV of the four muons. Here, the impact
parameter $\chi^{2}$ is defined as the difference in the vertex-fit $\chi^{2}$
of a given PV reconstructed with and without the track or candidate being
considered. In Run 2, two variables reflecting the isolation of the $B$
candidate from other tracks in the event are also included in the training:
the changes in the decay vertex $\chi^{2}$ when either one or two of the
closest additional tracks in the underlying event are added to the decay
vertex fit. Finally, the classifier response is transformed to give a uniform
distribution between 0 and 1 for simulated
${{B}^{0}_{s}}\\!\rightarrow{\mu^{+}\mu^{-}\mu^{+}\mu^{-}}$ candidates.

In the searches for
${{B}_{({s})}^{0}}\\!\rightarrow{\mu^{+}\mu^{-}\mu^{+}\mu^{-}}$ and
${{B}_{({s})}^{0}}\\!\rightarrow{{J\mskip-3.0mu/\mskip-2.0mu\psi}}\left({\mu^{+}\mu^{-}}\right){\mu^{+}\mu^{-}}$
decays the data are split into four intervals in the BDT classifier response
with boundaries 0.2, 0.3, 0.5, 0.6 and 1.0, and a simultaneous fit to all four
BDT intervals is used to search for the decays. Due to the low background
levels in the
${{B}_{({s})}^{0}}\\!\rightarrow{a}\left({\mu^{+}\mu^{-}}\right){a}\left({\mu^{+}\mu^{-}}\right)$
search region, the BDT response in this case is only required to be greater
than 0.1. Similarly, the selection for
${{B}^{0}_{s}}\\!\rightarrow{{J\mskip-3.0mu/\mskip-2.0mu\psi}}\left({\mu^{+}\mu^{-}}\right)\phi\left({\mu^{+}\mu^{-}}\right)$
decays requires the BDT classifier response to exceed 0.05.

##  4 Background

A range of background sources have the potential to contaminate the signal
search regions. The four-muon decays
${{B}^{0}_{s}}\\!\rightarrow{{J\mskip-3.0mu/\mskip-2.0mu\psi}}\left({\mu^{+}\mu^{-}}\right)\phi\left({\mu^{+}\mu^{-}}\right)$
and
${{B}^{0}_{s}}\\!\rightarrow\phi\left({\mu^{+}\mu^{-}}\right){\mu^{+}\mu^{-}}$,
which have branching fractions of $1.74\times 10^{-8}$ [19] and $2.3\times
10^{-10}$ [18], respectively, are reduced to negligible levels by the $\phi$
and ${J\mskip-3.0mu/\mskip-2.0mu\psi}$ vetoes. Potential background to the
${{B}_{({s})}^{0}}\\!\rightarrow{a}\left({\mu^{+}\mu^{-}}\right){a}\left({\mu^{+}\mu^{-}}\right)$
search comes from
${{B}^{0}_{s}}\\!\rightarrow\phi\left({\mu^{+}\mu^{-}}\right)\phi\left({\mu^{+}\mu^{-}}\right)$
decays, however the corresponding branching fraction of $1.5\times 10^{-12}$
is well below the sensitivity of this search.

A second category of background arises from decays of $b$-hadrons of the form
$H_{b}\rightarrow h^{+}h^{{}^{\prime}-}{\mu^{+}}{\mu^{-}}$ where two hadrons
($h^{+}$ and $h^{{}^{\prime}-}$) are misidentified as muons. Such background
includes
${{B}^{0}_{s}}\\!\rightarrow{{J\mskip-3.0mu/\mskip-2.0mu\psi}}\left({\mu^{+}\mu^{-}}\right)\phi\left({{K}^{+}}{{K}^{-}}\right)$,
${{B}^{0}_{s}}\\!\rightarrow{{J\mskip-3.0mu/\mskip-2.0mu\psi}}\left({\mu^{+}\mu^{-}}\right){{K}^{+}}{{K}^{-}}$,
${{B}^{0}_{s}}\\!\rightarrow{{J\mskip-3.0mu/\mskip-2.0mu\psi}}\left({\mu^{+}\mu^{-}}\right){{\pi}^{+}}{{\pi}^{-}}$,
${{B}^{0}}\\!\rightarrow{{J\mskip-3.0mu/\mskip-2.0mu\psi}}\left({\mu^{+}\mu^{-}}\right){{K}^{*0}}\left({{K}^{+}}{{\pi}^{-}}\right)$,
${{B}^{0}}\\!\rightarrow{{K}^{*0}}\left({{K}^{+}}{{\pi}^{-}}\right){\mu^{+}\mu^{-}}$,
${{B}^{0}}\\!\rightarrow\rho\left({{\pi}^{+}}{{\pi}^{-}}\right){\mu^{+}\mu^{-}}$,
${{\mathchar 28931\relax}^{0}_{b}}\\!\rightarrow\mathchar
28931\relax\left({p}{{\pi}^{-}}\right){\mu^{+}\mu^{-}}$ and ${{\mathchar
28931\relax}^{0}_{b}}\\!\rightarrow\mathchar
28931\relax\left({p}{{K}^{-}}\right){\mu^{+}\mu^{-}}$ decays. The typical
probability for a hadron to be misidentified as a muon by the LHCb detector is
around 1% [31], and this rate is further reduced by the PID requirements
described in Sec. 3. Studies performed using simulation demonstrate that all
these decays are reduced to negligible levels by the signal and the control
sample selections. As such, this analysis is essentially free of background
from other $b$-hadron decays, which are therefore not modelled in the fits for
the signal and control decays.

##  5 Fit procedure

The signal yields in the branching fraction measurements are normalised with
respect to the
${{B}^{0}_{s}}\\!\rightarrow{{J\mskip-3.0mu/\mskip-2.0mu\psi}}\left({\mu^{+}\mu^{-}}\right)\phi\left({\mu^{+}\mu^{-}}\right)$
control mode such that the branching fraction of a given signal decay is
calculated as

$\mathcal{B}_{\mathrm{sig}}=\alpha^{n}_{s/d}\times N^{n}_{\mathrm{sig}},$ (2)

where the normalisation factors for ${B}^{0}_{s}$ and ${B}^{0}$ modes,
$\alpha^{n}_{s}$ and $\alpha^{n}_{d}$, are defined as

$\displaystyle\alpha^{n}_{s}$
$\displaystyle\equiv\frac{\mathcal{B}_{\text{norm}}}{N_{\text{norm}}}\cdot\frac{\epsilon_{\text{norm}}}{\epsilon_{\text{sig}}^{n}},$
$\displaystyle\alpha^{n}_{d}$
$\displaystyle\equiv\frac{\mathcal{B}_{\text{norm}}}{N_{\text{norm}}}\cdot\frac{\epsilon_{\text{norm}}}{\epsilon_{\text{sig}}^{n}}\cdot\frac{f_{s}}{f_{d}}.$
(3)

The symbols sig and norm refer to the signal and normalisation decays,
respectively, $\mathcal{B}$ is the branching fraction, $N$ is the yield,
$\epsilon$ is the total efficiency to reconstruct and select a given decay
mode and $n$ indexes the BDT response intervals used in the simultaneous fits
for ${{B}_{({s})}^{0}}\\!\rightarrow{\mu^{+}\mu^{-}\mu^{+}\mu^{-}}$ and
${{B}_{({s})}^{0}}\\!\rightarrow{{J\mskip-3.0mu/\mskip-2.0mu\psi}}\left({\mu^{+}\mu^{-}}\right){\mu^{+}\mu^{-}}$
decays (in the case of
${{B}_{({s})}^{0}}\\!\rightarrow{a}\left({\mu^{+}\mu^{-}}\right){a}\left({\mu^{+}\mu^{-}}\right)$
there is only a single BDT response interval).

The branching fraction of the
${{B}^{0}_{s}}\\!\rightarrow{{J\mskip-3.0mu/\mskip-2.0mu\psi}}\left({\mu^{+}\mu^{-}}\right)\phi\left({\mu^{+}\mu^{-}}\right)$
normalisation mode is calculated from the product of the branching fractions
of the ${{B}^{0}_{s}}\\!\rightarrow{{J\mskip-3.0mu/\mskip-2.0mu\psi}}\phi$
[19], ${{J\mskip-3.0mu/\mskip-2.0mu\psi}}\\!\rightarrow{\mu^{+}\mu^{-}}$ and
$\phi\\!\rightarrow{\mu^{+}\mu^{-}}$ [18] decays, yielding
$\mathcal{B}_{\mathrm{norm}}=\left(1.74\pm 0.14\right)\times 10^{-8}$. The
ratio of fragmentation fractions is calculated as a weighted average of values
measured at centre-of-mass energies of 7, 8 and 13 TeV [19], giving
$f_{s}/f_{d}=0.250\pm 0.008$. The correlation in the measurements of
$\mathcal{B}\left({{B}^{0}_{s}}\\!\rightarrow{{J\mskip-3.0mu/\mskip-2.0mu\psi}}\phi\right)$
and $f_{s}/f_{d}$ is taken into account when calculating the uncertainty in
the normalisation term $\alpha_{d}^{i}$.

The efficiencies are calculated using simulated decays, to which weights are
applied in order to improve concordance with data. These weights are
calculated by comparing
${{B}^{0}_{s}}\\!\rightarrow{{J\mskip-3.0mu/\mskip-2.0mu\psi}}\left({\mu^{+}\mu^{-}}\right)\phi\left({{K}^{+}}{{K}^{-}}\right)$
decays in data and simulation, for which the branching fraction is three
orders of magnitude larger than for
${{B}^{0}_{s}}\\!\rightarrow{{J\mskip-3.0mu/\mskip-2.0mu\psi}}\left({\mu^{+}\mu^{-}}\right)\phi\left({\mu^{+}\mu^{-}}\right)$
and therefore allows for a much more precise determination of differences
between data and simulation. The trigger and offline selections for these
decays are similar to those used for
${{B}^{0}_{s}}\\!\rightarrow{{J\mskip-3.0mu/\mskip-2.0mu\psi}}\left({\mu^{+}\mu^{-}}\right)\phi\left({\mu^{+}\mu^{-}}\right)$
decays, but no PID requirement is applied on the kaons. The distributions of
the variables of interest for
${{B}^{0}_{s}}\\!\rightarrow{{J\mskip-3.0mu/\mskip-2.0mu\psi}}\left({\mu^{+}\mu^{-}}\right)\phi\left({{K}^{+}}{{K}^{-}}\right)$
decays are separated from background in data using the sPlot method [34],
using the $B$ candidate invariant mass as the discriminating variable.

A first set of weights, referred to as the generator weights, corrects the
$p_{\mathrm{T}}$ and pseudorapidity distributions of $B$ mesons as generated
by Pythia 8 [22], along with the multiplicity of the underlying events. A
second set of weights, referred to as the reconstruction weights, is used to
correct the vertex $\chi^{2}$ per degree of freedom and the impact parameter
$\chi^{2}$ of the $B$ mesons. The efficiencies for each mode in each data-
taking year are then calculated as

$\epsilon=\epsilon_{\text{presel}}\cdot\epsilon_{\text{BDT}}=\frac{\sum\limits_{\text{presel}}\omega_{\text{gen}}}{\sum\limits_{\text{gen}}\omega_{\text{gen}}}\cdot\frac{\sum\limits_{\text{BDT}}\omega_{\text{gen}}\omega_{\text{rec}}}{\sum\limits_{\text{presel}}\omega_{\text{gen}}\omega_{\text{rec}}}\,,$
(4)

where $\sum\limits_{\text{gen}}$, $\sum\limits_{\text{presel}}$ and
$\sum\limits_{\text{BDT}}$ refer to sums over the generated simulation sample
before any reconstruction or selection, the offline reconstructed simulation
sample passing the full selection excluding the BDT requirements, and the
offline reconstructed simulation sample passing the full selection including
the BDT requirements, respectively. Meanwhile, $\omega_{\text{gen}}$ and
$\omega_{\text{rec}}$ refer to the generator and reconstruction weights. Note
that the reconstruction weights are only used in the calculation of the
efficiencies of the BDT requirements, since requirements on the $B$
candidate’s vertex $\chi^{2}$ per degree of freedom and the impact parameter
$\chi^{2}$ imposed in prior stages of the selection are highly efficient for
the signal modes. The individual efficiencies for each year are then combined
in a weighted average according to the integrated luminosity recorded and the
$b\overline{b}$ production cross-section in that year [35, 36].

The invariant mass distribution of
${{B}^{0}_{s}}\\!\rightarrow{{J\mskip-3.0mu/\mskip-2.0mu\psi}}\left({\mu^{+}\mu^{-}}\right)\phi\left({\mu^{+}\mu^{-}}\right)$
candidates in the range
$5100<m({\mu^{+}\mu^{-}}{\mu^{+}\mu^{-}})<6000\text{\,Me\kern-1.00006ptV\\!/}c^{2}$
is shown in Fig. 3. The yield of the normalisation mode is determined using an
extended unbinned maximum-likelihood fit to this invariant mass spectrum. The
${{B}^{0}_{s}}\\!\rightarrow{{J\mskip-3.0mu/\mskip-2.0mu\psi}}\left({\mu^{+}\mu^{-}}\right)\phi\left({\mu^{+}\mu^{-}}\right)$
signal is modelled using the sum of two Crystal Ball functions [37] with
common mean, referred to as a double Crystal Ball (DCB) function. The
parameters of this function are determined from a fit to simulated
${{B}^{0}_{s}}\\!\rightarrow{{J\mskip-3.0mu/\mskip-2.0mu\psi}}\left({\mu^{+}\mu^{-}}\right)\phi\left({\mu^{+}\mu^{-}}\right)$
decays, in which the values are fixed in the fit to data, with the exception
of the mean and widths of the DCB function, which are allowed to vary freely
in order to account for differences in resolution and mass scale between data
and simulation. The resulting ratio of the mass resolution in data over
simulation is found to be $1.11\pm 0.09$. The combinatorial background is
modelled with an exponential function. Background from the six signal decays
in the control mode mass distribution is negligible due to the low branching
fractions of the signal decays and the requirements on the
${J\mskip-3.0mu/\mskip-2.0mu\psi}$ and $\phi$ masses imposed in the control
mode selection. The yield of
${{B}^{0}_{s}}\\!\rightarrow{{J\mskip-3.0mu/\mskip-2.0mu\psi}}\left({\mu^{+}\mu^{-}}\right)\phi\left({\mu^{+}\mu^{-}}\right)$
decays is found to be $218\pm 16$, where the uncertainty is statistical only.
The invariant mass spectrum and the fit projection are shown in Fig. 3.

Figure 3: Distribution of the ${\mu^{+}\mu^{-}}{\mu^{+}\mu^{-}}$ invariant
mass of candidates passing the
${{B}^{0}_{s}}\\!\rightarrow{{J\mskip-3.0mu/\mskip-2.0mu\psi}}\left({\mu^{+}\mu^{-}}\right)\phi\left({\mu^{+}\mu^{-}}\right)$
selection, with the fit model used to determine the control sample yield
overlaid.

As described in Sec. 3, the data used to search for
${{B}_{({s})}^{0}}\\!\rightarrow{\mu^{+}\mu^{-}\mu^{+}\mu^{-}}$ and
${{B}_{({s})}^{0}}\\!\rightarrow{{J\mskip-3.0mu/\mskip-2.0mu\psi}}\left({\mu^{+}\mu^{-}}\right){\mu^{+}\mu^{-}}$
decays are split into four intervals in the BDT classifier response, while for
the
${{B}_{({s})}^{0}}\\!\rightarrow{a}\left({\mu^{+}\mu^{-}}\right){a}\left({\mu^{+}\mu^{-}}\right)$
data sample a single BDT interval is defined. The normalisation factors
$\alpha_{s}^{i}$ and $\alpha_{d}^{i}$ are calculated using Eq. 2 and the
inputs described above, with their values Gaussian constrained to their
central values according to their statistical uncertainties. Additional
Gaussian constraints are used to include the effects of the systematic
uncertainties in the inputs described in Sec. 6.

For each of the
${{B}_{({s})}^{0}}\\!\rightarrow{\mu^{+}\mu^{-}\mu^{+}\mu^{-}}$ and
${{B}_{({s})}^{0}}\\!\rightarrow{{J\mskip-3.0mu/\mskip-2.0mu\psi}}\left({\mu^{+}\mu^{-}}\right){\mu^{+}\mu^{-}}$
signal modes, an unbinned extended maximum-likelihood fit is performed to the
$B$ candidate invariant mass distribution in the range
$4900<m({\mu^{+}\mu^{-}}{\mu^{+}\mu^{-}})<6000\text{\,Me\kern-1.00006ptV\\!/}c^{2}$,
simultaneously in all four BDT intervals, in order to search for the decays
and measure their branching fractions. In the case of the
${{B}_{({s})}^{0}}\\!\rightarrow{a}\left({\mu^{+}\mu^{-}}\right){a}\left({\mu^{+}\mu^{-}}\right)$
modes, fits are performed to the invariant mass distribution in a single
region of the BDT response. All signal mass shapes are modelled using DCB
functions, in which the parameters are fixed to values obtained from
simulation, with the exception of the mean and widths, which are offset and
scaled respectively according to the results of the fit for the control mode.
Due to their overlap, the fits to the ${B}^{0}_{s}$ and ${B}^{0}$ decays are
performed separately in each case, i.e. the ${B}^{0}_{s}$ component is set to
zero when performing the fit to the ${B}^{0}$ component and vice versa. The
yields of each signal mode are not fitted directly, rather they are expressed
in terms of the normalisation factors and the decay branching fractions (see
Eq. 2), which are allowed to vary freely. The branching fraction is
constrained to be greater than or equal to zero in each fit. The combinatorial
background is modelled using exponential functions, in which the slope
parameters and yields are allowed to vary freely independently in each BDT
interval.

##  6 Systematic uncertainties

Due to the similarity between the signal and control decay modes, many sources
of systematic uncertainty cancel in the ratio of efficiencies. Nonetheless, a
number of sources of systematic uncertainty have the potential to affect the
measured branching fractions. Two sources of uncertainty arise from the models
used to generate the simulated signal decays. In particular, since no
theoretical prediction for the decay dynamics of
${{B}_{({s})}^{0}}\\!\rightarrow{\mu^{+}\mu^{-}\mu^{+}\mu^{-}}$ decays is
currently available, decays are generated using a phase-space model. This may
lead to a bias in the efficiency calculation if the true kinematic
distributions of these decays differs from the one used in the simulation. To
estimate the representative scale of the bias arising from such an effect, the
dimuon invariant mass squared ($q^{2}_{ij}$) distributions for
${{B}_{({s})}^{0}}\\!\rightarrow{\mu^{+}\mu^{-}\mu^{+}\mu^{-}}$ decays are
weighted to be uniform in the $(q^{2}_{11},q^{2}_{22})$ plane, and the
efficiencies are evaluated. A similar effect is evaluated for
${{B}_{({s})}^{0}}\\!\rightarrow{{J\mskip-3.0mu/\mskip-2.0mu\psi}}\left({\mu^{+}\mu^{-}}\right){\mu^{+}\mu^{-}}$
decays, where the invariant mass squared of the two muons not originating from
the ${J\mskip-3.0mu/\mskip-2.0mu\psi}$ meson is weighted to give a uniform
distribution. This leads to shifts relative to the nominal efficiencies of up
to about 20% depending on the decay mode and BDT interval, which are taken as
systematic uncertainties. This represents the largest source of systematic
uncertainty. Similarly, the effective lifetimes of the ${B}^{0}_{s}$ meson
decays, ${{B}^{0}_{s}}\\!\rightarrow{\mu^{+}\mu^{-}\mu^{+}\mu^{-}}$,
${{B}^{0}_{s}}\\!\rightarrow{a}\left({\mu^{+}\mu^{-}}\right){a}\left({\mu^{+}\mu^{-}}\right)$
and
${{B}^{0}_{s}}\\!\rightarrow{{J\mskip-3.0mu/\mskip-2.0mu\psi}}\left({\mu^{+}\mu^{-}}\right){\mu^{+}\mu^{-}}$,
depend on the CP-even and CP-odd contributions in the decay amplitudes, which
are a priori unknown. These decays are generated assuming the mean
${B}^{0}_{s}$ lifetime. If their true effective lifetimes differ significantly
from this value, then the efficiency estimates will be biased. The maximum
size of such an effect is estimated by weighting simulated decay modes to have
their effective lifetimes equal to the lifetimes of the heavy and light
${B}^{0}_{s}$ mass eigenstates. The largest shift with respect to the nominal
efficiency is again taken as a source of systematic uncertainty, with relative
shifts in the efficiencies at around the 5% level.

Further sources of uncertainty arise due to data-simulation differences. The
effect of mismodelling of the PID response in simulation is evaluated using
calibration samples of muons from data to compare data-driven efficiencies
with those obtained directly from simulation, with the differences in the
ratios of efficiencies of 1–2% taken as systematic uncertainties. The effect
of the difference in mass resolution between simulation and data on the
selection efficiency is estimated to be below 1% using the results of the fits
to the control mode mass distribution. Finally, the efficiencies are evaluated
without the application of the generator and reconstruction weights described
in Sec. 5 in order to get an upper bound on the likely effect of mismodelling
in simulation. Since the effect on the efficiency ratio is found to be small,
up to around 10% depending on the decay and BDT interval, this shift is
conservatively taken as a source of systematic uncertainty in the efficiency.
All these systematic uncertainties are included in the branching fraction fit
by applying Gaussian constraints to the efficiencies or efficiency ratios.

##  7 Results

The invariant mass distributions of ${B}^{0}_{s}$ and ${B}^{0}$ candidates
passing the ${{B}_{({s})}^{0}}\\!\rightarrow{\mu^{+}\mu^{-}\mu^{+}\mu^{-}}$,
${{B}_{({s})}^{0}}\\!\rightarrow{a}\left({\mu^{+}\mu^{-}}\right){a}\left({\mu^{+}\mu^{-}}\right)$
and
${{B}_{({s})}^{0}}\\!\rightarrow{{J\mskip-3.0mu/\mskip-2.0mu\psi}}\left({\mu^{+}\mu^{-}}\right){\mu^{+}\mu^{-}}$
selections in the most sensitive BDT intervals are shown in Figs. 4, 5 and 6,
respectively. The projections of the fits used to search for the signal decays
and measure their branching fractions are overlaid. The distributions and fit
projections in the remaining BDT intervals are shown in Appendix A. No
evidence for any of the six signal decay modes is found, with the most
significant excesses found in the
${{B}^{0}_{s}}\\!\rightarrow{{J\mskip-3.0mu/\mskip-2.0mu\psi}}\left({\mu^{+}\mu^{-}}\right){\mu^{+}\mu^{-}}$
and
${{B}^{0}}\\!\rightarrow{{J\mskip-3.0mu/\mskip-2.0mu\psi}}\left({\mu^{+}\mu^{-}}\right){\mu^{+}\mu^{-}}$
searches, amounting to two standard deviations in both cases, calculated using
Wilks’ theorem [38]. Limits are set on the branching fractions using the
$\text{CL}_{\text{s}}$ method [39] as implemented in the GammaCombo package
[40, 41] using a one-sided test statistic, with 80 scan points and 2000
pseudoexperiments per scan point. The limits at 95% confidence level are

$\displaystyle\mathcal{B}\left({{B}^{0}_{s}}\\!\rightarrow{\mu^{+}\mu^{-}\mu^{+}\mu^{-}}\right)$
$\displaystyle<\,$ $\displaystyle 8.6\times 10^{-10}\,,$
$\displaystyle\mathcal{B}\left({{B}^{0}}\\!\rightarrow{\mu^{+}\mu^{-}\mu^{+}\mu^{-}}\right)$
$\displaystyle<\,$ $\displaystyle 1.8\times 10^{-10}\,,$
$\displaystyle\mathcal{B}\left({{B}^{0}_{s}}\\!\rightarrow{a}\left({\mu^{+}\mu^{-}}\right){a}\left({\mu^{+}\mu^{-}}\right)\right)$
$\displaystyle<\,$ $\displaystyle 5.8\times 10^{-10}\,,$
$\displaystyle\mathcal{B}\left({{B}^{0}}\\!\rightarrow{a}\left({\mu^{+}\mu^{-}}\right){a}\left({\mu^{+}\mu^{-}}\right)\right)$
$\displaystyle<\,$ $\displaystyle 2.3\times 10^{-10}\,,$
$\displaystyle\mathcal{B}\left({{B}^{0}_{s}}\\!\rightarrow{{J\mskip-3.0mu/\mskip-2.0mu\psi}}\left({\mu^{+}\mu^{-}}\right){\mu^{+}\mu^{-}}\right)$
$\displaystyle<\,$ $\displaystyle 2.6\times 10^{-9}\,,$
$\displaystyle\mathcal{B}\left({{B}^{0}}\\!\rightarrow{{J\mskip-3.0mu/\mskip-2.0mu\psi}}\left({\mu^{+}\mu^{-}}\right){\mu^{+}\mu^{-}}\right)$
$\displaystyle<\,$ $\displaystyle 1.0\times 10^{-9}\,.$

The corresponding $\text{CL}_{\text{s}}$ scans for all six decays are shown in
Fig. 7. Note that in the case of
${{B}_{({s})}^{0}}\\!\rightarrow{a}\left({\mu^{+}\mu^{-}}\right){a}\left({\mu^{+}\mu^{-}}\right)$
decays, these limits are evaluated assuming a promptly decaying intermediate
scalar with a mass of 1$\text{\,Ge\kern-1.00006ptV\\!/}c^{2}$. The limits
quoted for
${{B}_{({s})}^{0}}\\!\rightarrow{{J\mskip-3.0mu/\mskip-2.0mu\psi}}\left({\mu^{+}\mu^{-}}\right){\mu^{+}\mu^{-}}$
decays include the
${{J\mskip-3.0mu/\mskip-2.0mu\psi}}\\!\rightarrow{\mu^{+}\mu^{-}}$ branching
fraction. These results constitute the most stringent limits on each of the
six ${{B}_{({s})}^{0}}\\!\rightarrow{\mu^{+}\mu^{-}\mu^{+}\mu^{-}}$ decays to
date and supersede previous results by the LHCb experiment [16].

Figure 4: Distribution of the ${\mu^{+}\mu^{-}}{\mu^{+}\mu^{-}}$ invariant
mass of candidates passing the
${{B}_{({s})}^{0}}\\!\rightarrow{\mu^{+}\mu^{-}\mu^{+}\mu^{-}}$ selection in
the most sensitive BDT interval, with the fit models used to determine the
branching fractions of (left)
${{B}^{0}_{s}}\\!\rightarrow{\mu^{+}\mu^{-}\mu^{+}\mu^{-}}$ and (right)
${{B}^{0}}\\!\rightarrow{\mu^{+}\mu^{-}\mu^{+}\mu^{-}}$ decays overlaid.

Figure 5: Distribution of the ${\mu^{+}\mu^{-}}{\mu^{+}\mu^{-}}$ invariant
mass of candidates passing the
${{B}_{({s})}^{0}}\\!\rightarrow{a}\left({\mu^{+}\mu^{-}}\right){a}\left({\mu^{+}\mu^{-}}\right)$
selection, with the fit models used to determine the branching fractions of
(left)
${{B}^{0}_{s}}\\!\rightarrow{a}\left({\mu^{+}\mu^{-}}\right){a}\left({\mu^{+}\mu^{-}}\right)$
and (right)
${{B}^{0}}\\!\rightarrow{a}\left({\mu^{+}\mu^{-}}\right){a}\left({\mu^{+}\mu^{-}}\right)$
decays overlaid.

Figure 6: Distribution of the ${\mu^{+}\mu^{-}}{\mu^{+}\mu^{-}}$ invariant
mass of candidates passing the
${{B}_{({s})}^{0}}\\!\rightarrow{{J\mskip-3.0mu/\mskip-2.0mu\psi}}\left({\mu^{+}\mu^{-}}\right){\mu^{+}\mu^{-}}$
selection in the most sensitive BDT interval, with the fit models used to
determine the branching fractions of (left)
${{B}^{0}_{s}}\\!\rightarrow{{J\mskip-3.0mu/\mskip-2.0mu\psi}}\left({\mu^{+}\mu^{-}}\right){\mu^{+}\mu^{-}}$
and (right)
${{B}^{0}}\\!\rightarrow{{J\mskip-3.0mu/\mskip-2.0mu\psi}}\left({\mu^{+}\mu^{-}}\right){\mu^{+}\mu^{-}}$
decays overlaid.

Figure 7: Confidence levels ($\mathrm{CL_{s}}$) obtained from the
$\text{CL}_{\text{s}}$ method to set limits on the branching fractions of (top
left) ${{B}^{0}_{s}}\\!\rightarrow{\mu^{+}\mu^{-}\mu^{+}\mu^{-}}$, (top right)
${{B}^{0}}\\!\rightarrow{\mu^{+}\mu^{-}\mu^{+}\mu^{-}}$, (middle left)
${{B}^{0}_{s}}\\!\rightarrow{a}\left({\mu^{+}\mu^{-}}\right){a}\left({\mu^{+}\mu^{-}}\right)$,
(middle right)
${{B}^{0}}\\!\rightarrow{a}\left({\mu^{+}\mu^{-}}\right){a}\left({\mu^{+}\mu^{-}}\right)$,
(bottom left)
${{B}^{0}_{s}}\\!\rightarrow{{J\mskip-3.0mu/\mskip-2.0mu\psi}}\left({\mu^{+}\mu^{-}}\right){\mu^{+}\mu^{-}}$
and (bottom right)
${{B}^{0}}\\!\rightarrow{{J\mskip-3.0mu/\mskip-2.0mu\psi}}\left({\mu^{+}\mu^{-}}\right){\mu^{+}\mu^{-}}$
decays.

## Acknowledgements

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); MICINN (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 NERSC (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 ARC and ARDC (Australia); AvH Foundation
(Germany); EPLANET, Marie Skłodowska-Curie Actions and ERC (European Union);
A*MIDEX, ANR, IPhU and Labex P2IO, and Région Auvergne-Rhône-Alpes (France);
Key Research Program of Frontier Sciences of CAS, CAS PIFI, CAS CCEPP,
Fundamental Research Funds for the Central Universities, and Sci. & Tech.
Program of Guangzhou (China);

RFBR, RSF and Yandex LLC (Russia); GVA, XuntaGal and GENCAT (Spain); the
Leverhulme Trust, the Royal Society and UKRI (United Kingdom).

##  Appendix A Additional figures

This appendix contains the distributions of ${\mu^{+}\mu^{-}}{\mu^{+}\mu^{-}}$
invariant mass for candidates passing the
${{B}_{({s})}^{0}}\\!\rightarrow{\mu^{+}\mu^{-}\mu^{+}\mu^{-}}$ and
${{B}_{({s})}^{0}}\\!\rightarrow{{J\mskip-3.0mu/\mskip-2.0mu\psi}}\left({\mu^{+}\mu^{-}}\right){\mu^{+}\mu^{-}}$
selections, in the lower three intervals of the BDT classifier response, with
the results of the maximum-likelihood fits overlaid.

Figure 8: Distribution of the ${\mu^{+}\mu^{-}}{\mu^{+}\mu^{-}}$ invariant
mass of candidates passing the
${{B}_{({s})}^{0}}\\!\rightarrow{\mu^{+}\mu^{-}\mu^{+}\mu^{-}}$ selection in
(top left) the lowest BDT interval, (top right) the second lowest BDT interval
and (bottom) the second highest BDT interval, with the fit models used to
determine the branching fraction of
${{B}^{0}_{s}}\\!\rightarrow{\mu^{+}\mu^{-}\mu^{+}\mu^{-}}$ overlaid.

Figure 9: Distribution of the ${\mu^{+}\mu^{-}}{\mu^{+}\mu^{-}}$ invariant
mass of candidates passing the
${{B}_{({s})}^{0}}\\!\rightarrow{\mu^{+}\mu^{-}\mu^{+}\mu^{-}}$ selection in
(top left) the lowest BDT interval, (top right) the second lowest BDT interval
and (bottom) the second highest BDT interval, with the fit models used to
determine the branching fraction of
${{B}^{0}}\\!\rightarrow{\mu^{+}\mu^{-}\mu^{+}\mu^{-}}$ overlaid.

Figure 10: Distribution of the ${\mu^{+}\mu^{-}}{\mu^{+}\mu^{-}}$ invariant
mass of candidates passing the
${{B}_{({s})}^{0}}\\!\rightarrow{{J\mskip-3.0mu/\mskip-2.0mu\psi}}\left({\mu^{+}\mu^{-}}\right){\mu^{+}\mu^{-}}$
selection in (top left) the lowest BDT interval, (top right) the second lowest
BDT interval and (bottom) the second highest BDT interval, with the fit models
used to determine the branching fraction of
${{B}^{0}_{s}}\\!\rightarrow{{J\mskip-3.0mu/\mskip-2.0mu\psi}}\left({\mu^{+}\mu^{-}}\right){\mu^{+}\mu^{-}}$
overlaid.

Figure 11: Distribution of the ${\mu^{+}\mu^{-}}{\mu^{+}\mu^{-}}$ invariant
mass of candidates passing the
${{B}_{({s})}^{0}}\\!\rightarrow{{J\mskip-3.0mu/\mskip-2.0mu\psi}}\left({\mu^{+}\mu^{-}}\right){\mu^{+}\mu^{-}}$
selection in (top left) the lowest BDT interval, (top right) the second lowest
BDT interval and (bottom) the second highest BDT interval, with the fit models
used to determine the branching fraction of
${{B}^{0}}\\!\rightarrow{{J\mskip-3.0mu/\mskip-2.0mu\psi}}\left({\mu^{+}\mu^{-}}\right){\mu^{+}\mu^{-}}$
overlaid.

## References

  * [1] A. V. Danilina and N. V. Nikitin, _Four-leptonic decays of charged and neutral $B$ mesons within the Standard Model_, Phys. Atom. Nucl. 81 (2018) 347
  * [2] S. V. Demidov and D. S. Gorbunov, _Flavor violating processes with sgoldstino pair production_ , Phys. Rev. D85 (2012) 077701, arXiv:1112.5230
  * [3] M. Bauer, M. Neubert, and A. Thamm, _LHC as an axion factory: probing an axion explanation for $(g-2)_{\mu}$ with exotic Higgs decays_, Phys. Rev. Lett. 119 (2017) 031802, arXiv:1704.08207
  * [4] J. Liu, C. E. M. Wagner, and X.-P. Wang, _A light complex scalar for the electron and muon anomalous magnetic moments_ , JHEP 03 (2019) 008, arXiv:1810.11028
  * [5] T. Aoyama et al., _The anomalous magnetic moment of the muon in the Standard Model_ , Phys. Rept. 887 (2020) 1, arXiv:2006.04822
  * [6] Muon g–2 collaboration, B. Abi et al., _Measurement of the positive muon anomalous magnetic moment to 0.46 ppm_ , Phys. Rev. Lett. 126 (2021) 141801, arXiv:2104.03281
  * [7] Muon g–2 collaboration, G. W. Bennett et al., _Final report of the muon E821 anomalous magnetic moment measurement at BNL_ , Phys. Rev. D73 (2006) 072003, arXiv:hep-ex/0602035
  * [8] LHCb collaboration, R. Aaij et al., _Tests of lepton universality using $B^{0}\rightarrow K^{0}_{S}\ell^{+}\ell^{-}$ and $B^{+}\rightarrow K^{*+}\ell^{+}\ell^{-}$ decays_, arXiv:2110.09501, submitted to PRL
  * [9] LHCb collaboration, R. Aaij et al., _Test of lepton universality in beauty-quark decays_ , arXiv:2103.11769, to appear in Nature Physics
  * [10] LHCb collaboration, R. Aaij et al., _Test of lepton universality using ${{\mathchar 28931\relax}^{0}_{b}}\\!\rightarrow p{{K}^{-}}\ell^{+}\ell^{-}$ decays_, JHEP 05 (2020) 040, arXiv:1912.08139
  * [11] LHCb collaboration, R. Aaij et al., _Test of lepton universality with ${{B}^{0}}\\!\rightarrow{{K}^{*0}}\ell^{+}\ell^{-}$ decays_, JHEP 08 (2017) 055, arXiv:1705.05802
  * [12] LHCb collaboration, R. Aaij et al., _Angular analysis of the rare decay ${{B}^{0}_{s}}\rightarrow\phi{\mu^{+}}{\mu^{-}}$_, JHEP 11 (2021) 043, arXiv:2107.13428
  * [13] LHCb collaboration, R. Aaij et al., _Angular analysis of the $B^{+}\rightarrow K^{\ast+}\mu^{+}\mu^{-}$ decay_, Phys. Rev. Lett. 126 (2021) 161802, arXiv:2012.13241
  * [14] LHCb collaboration, R. Aaij et al., _Measurement of $C\\!P$-averaged observables in the ${{B}^{0}}\\!\rightarrow{{K}^{*0}}{\mu^{+}\mu^{-}}$ decay_, Phys. Rev. Lett. 125 (2020) 011802, arXiv:2003.04831
  * [15] M. Chala, U. Egede, and M. Spannowsky, _Searching new physics in rare $B$-meson decays into multiple muons_, Eur. Phys. J. C79 (2019) 431, arXiv:1902.10156
  * [16] LHCb collaboration, R. Aaij et al., _Search for decays of neutral beauty mesons into four muons_ , JHEP 03 (2017) 001, arXiv:1611.07704
  * [17] D. H. Evans, B. Grinstein, and D. R. Nolte, _Short distance analysis of ${{\kern 1.79993pt\overline{\kern-1.79993ptB}}{}^{0}}\rightarrow{{J\mskip-3.0mu/\mskip-2.0mu\psi}}{e^{+}}{e^{-}}$, ${{\kern 1.79993pt\overline{\kern-1.79993ptB}}{}^{0}}\rightarrow{\eta_{c}}{e^{+}}{e^{-}}$, ${{\kern 1.79993pt\overline{\kern-1.79993ptB}}{}^{0}}\rightarrow{{D}^{*0}}{e^{+}}{e^{-}}$ and ${{\kern 1.79993pt\overline{\kern-1.79993ptB}}{}^{0}}\rightarrow{{D}^{0}}{e^{+}}{e^{-}}$_, Nucl. Phys. B577 (2000) 240, arXiv:hep-ph/9906528
  * [18] Particle Data Group, P. A. Zyla et al., _Review of particle physics_ , Prog. Theor. Exp. Phys. 2020 (2020) 083C01
  * [19] LHCb collaboration, R. Aaij et al., _Precise measurement of the $f_{s}/f_{d}$ ratio of fragmentation fractions and of $B^{0}_{s}$ decay branching fractions_, Phys. Rev. D104 (2021) 032005, arXiv:2103.06810
  * [20] LHCb collaboration, A. A. Alves Jr. et al., _The LHCb detector at the LHC_, JINST 3 (2008) S08005
  * [21] LHCb collaboration, R. Aaij et al., _LHCb detector performance_ , Int. J. Mod. Phys. A30 (2015) 1530022, arXiv:1412.6352
  * [22] T. Sjöstrand, S. Mrenna, and P. Skands, _A brief introduction to PYTHIA 8.1_ , Comput. Phys. Commun. 178 (2008) 852, arXiv:0710.3820
  * [23] T. Sjöstrand, S. Mrenna, and P. Skands, _PYTHIA 6.4 physics and manual_ , JHEP 05 (2006) 026, arXiv:hep-ph/0603175
  * [24] I. Belyaev et al., _Handling of the generation of primary events in Gauss, the LHCb simulation framework_ , J. Phys. Conf. Ser. 331 (2011) 032047
  * [25] D. J. Lange, _The EvtGen particle decay simulation package_ , Nucl. Instrum. Meth. A462 (2001) 152
  * [26] N. Davidson, T. Przedzinski, and Z. Was, _PHOTOS interface in C++: Technical and physics documentation_ , Comp. Phys. Comm. 199 (2016) 86, arXiv:1011.0937
  * [27] Geant4 collaboration, J. Allison et al., _Geant4 developments and applications_ , IEEE Trans. Nucl. Sci. 53 (2006) 270
  * [28] Geant4 collaboration, S. Agostinelli et al., _Geant4: A simulation toolkit_ , Nucl. Instrum. Meth. A506 (2003) 250
  * [29] M. Clemencic et al., _The LHCb simulation application, Gauss: Design, evolution and experience_, J. Phys. Conf. Ser. 331 (2011) 032023
  * [30] A. Ali, P. Ball, L. T. Handoko, and G. Hiller, _A comparative study of the decays $B\rightarrow$ ($K$, $K^{*)}\ell^{+}\ell^{-}$ in standard model and supersymmetric theories_, Phys. Rev. D61 (2000) 074024, arXiv:hep-ph/9910221
  * [31] F. Archilli et al., _Performance of the muon identification at LHCb_ , JINST 8 (2013) P10020, arXiv:1306.0249
  * [32] H. Voss, A. Hoecker, J. Stelzer, and F. Tegenfeldt, _TMVA - Toolkit for Multivariate Data Analysis with ROOT_ , PoS ACAT (2007) 040
  * [33] A. Blum, A. Kalai, and J. Langford, _Beating the hold-out: bounds for k-fold and progressive cross-validation_ , in Proceedings of the twelfth annual conference on computational learning theory,  COLT ’99, (New York, NY, USA), 203, ACM, 1999
  * [34] M. Pivk and F. R. Le Diberder, _sPlot: A statistical tool to unfold data distributions_ , Nucl. Instrum. Meth. A555 (2005) 356, arXiv:physics/0402083
  * [35] LHCb collaboration, R. Aaij et al., _Measurement of the $b$-quark production cross-section in 7 and 13 TeV ${p}{p}$ collisions_, Phys. Rev. Lett. 118 (2017) 052002, Erratum ibid. 119 (2017) 169901, arXiv:1612.05140
  * [36] LHCb collaboration, R. Aaij et al., _Measurement of differential ${b}{\overline{b}}$ and ${c}{\overline{c}}$ dijet cross-sections in the forward region of $pp$ collisions at $\sqrt{s}=13$ TeV_, JHEP 02 (2021) 023, arXiv:2010.09437
  * [37] T. Skwarnicki, A study of the radiative cascade transitions between the Upsilon-prime and Upsilon resonances, PhD thesis, Institute of Nuclear Physics, Krakow, 1986, DESY-F31-86-02
  * [38] S. S. Wilks, _The large-sample distribution of the likelihood ratio for testing composite hypotheses_ , Ann. Math. Stat. 9 (1938) 60
  * [39] A. L. Read, _Presentation of search results: The CL s technique_, J. Phys. G28 (2002) 2693
  * [40] M. Kenzie et al., _GammaCombo: A statistical analysis framework for combining measurements, fitting datasets and producing confidence intervals_ ,  doi: 10.5281/zenodo.3371421
  * [41] LHCb collaboration, R. Aaij et al., _Measurement of the CKM angle $\gamma$ from a combination of LHCb results_, JHEP 12 (2016) 087, arXiv:1611.03076

LHCb collaboration

R. Aaij32, A.S.W. Abdelmotteleb56, C. Abellán Beteta50, F. Abudinén56, T.
Ackernley60, B. Adeva46, M. Adinolfi54, H. Afsharnia9, C. Agapopoulou13, C.A.
Aidala87, S. Aiola25, Z. Ajaltouni9, S. Akar65, J. Albrecht15, F. Alessio48,
M. Alexander59, A. Alfonso Albero45, Z. Aliouche62, G. Alkhazov38, P. Alvarez
Cartelle55, S. Amato2, J.L. Amey54, Y. Amhis11, L. An48, L. Anderlini22, M.
Andersson50, A. Andreianov38, M. Andreotti21, F. Archilli17, A. Artamonov44,
M. Artuso68, K. Arzymatov42, E. Aslanides10, M. Atzeni50, B. Audurier12, S.
Bachmann17, M. Bachmayer49, J.J. Back56, P. Baladron Rodriguez46, V.
Balagura12, W. Baldini21, J. Baptista Leite1, M. Barbetti22,h, R.J. Barlow62,
S. Barsuk11, W. Barter61, M. Bartolini55, F. Baryshnikov83, J.M. Basels14, S.
Bashir34, G. Bassi29, B. Batsukh68, A. Battig15, A. Bay49, A. Beck56, M.
Becker15, F. Bedeschi29, I. Bediaga1, A. Beiter68, V. Belavin42, S. Belin27,
V. Bellee50, K. Belous44, I. Belov40, I. Belyaev41, G. Bencivenni23, E. Ben-
Haim13, A. Berezhnoy40, R. Bernet50, D. Berninghoff17, H.C. Bernstein68, C.
Bertella62, A. Bertolin28, C. Betancourt50, F. Betti48, Ia. Bezshyiko50, S.
Bhasin54, J. Bhom35, L. Bian73, M.S. Bieker15, N.V. Biesuz21, S. Bifani53, P.
Billoir13, A. Biolchini32, M. Birch61, F.C.R. Bishop55, A. Bitadze62, A.
Bizzeti22,l, M. Bjørn63, M.P. Blago48, T. Blake56, F. Blanc49, S. Blusk68, D.
Bobulska59, J.A. Boelhauve15, O. Boente Garcia46, T. Boettcher65, A.
Boldyrev82, A. Bondar43, N. Bondar38,48, S. Borghi62, M. Borisyak42, M.
Borsato17, J.T. Borsuk35, S.A. Bouchiba49, T.J.V. Bowcock60,48, A. Boyer48, C.
Bozzi21, M.J. Bradley61, S. Braun66, A. Brea Rodriguez46, J. Brodzicka35, A.
Brossa Gonzalo56, D. Brundu27, A. Buonaura50, L. Buonincontri28, A.T. Burke62,
C. Burr48, A. Bursche72, A. Butkevich39, J.S. Butter32, J. Buytaert48, W.
Byczynski48, S. Cadeddu27, H. Cai73, R. Calabrese21,g, L. Calefice15,13, S.
Cali23, R. Calladine53, M. Calvi26,k, M. Calvo Gomez85, P. Camargo
Magalhaes54, P. Campana23, A.F. Campoverde Quezada6, S. Capelli26,k, L.
Capriotti20,e, A. Carbone20,e, G. Carboni31,q, R. Cardinale24,i, A. Cardini27,
I. Carli4, P. Carniti26,k, L. Carus14, K. Carvalho Akiba32, A. Casais Vidal46,
R. Caspary17, G. Casse60, M. Cattaneo48, G. Cavallero48, S. Celani49, J.
Cerasoli10, D. Cervenkov63, A.J. Chadwick60, M.G. Chapman54, M. Charles13, Ph.
Charpentier48, G. Chatzikonstantinidis53, C.A. Chavez Barajas60, M.
Chefdeville8, C. Chen3, S. Chen4, A. Chernov35, V. Chobanova46, S. Cholak49,
M. Chrzaszcz35, A. Chubykin38, V. Chulikov38, P. Ciambrone23, M.F. Cicala56,
X. Cid Vidal46, G. Ciezarek48, P.E.L. Clarke58, M. Clemencic48, H.V. Cliff55,
J. Closier48, J.L. Cobbledick62, V. Coco48, J.A.B. Coelho11, J. Cogan10, E.
Cogneras9, L. Cojocariu37, P. Collins48, T. Colombo48, L. Congedo19,d, A.
Contu27, N. Cooke53, G. Coombs59, I. Corredoira 46, G. Corti48, C.M. Costa
Sobral56, B. Couturier48, D.C. Craik64, J. Crkovská67, M. Cruz Torres1, R.
Currie58, C.L. Da Silva67, S. Dadabaev83, L. Dai71, E. Dall’Occo15, J.
Dalseno46, C. D’Ambrosio48, A. Danilina41, P. d’Argent48, A. Dashkina83, J.E.
Davies62, A. Davis62, O. De Aguiar Francisco62, K. De Bruyn79, S. De Capua62,
M. De Cian49, E. De Lucia23, J.M. De Miranda1, L. De Paula2, M. De Serio19,d,
D. De Simone50, P. De Simone23, F. De Vellis15, J.A. de Vries80, C.T. Dean67,
F. Debernardis19,d, D. Decamp8, V. Dedu10, L. Del Buono13, B. Delaney55, H.-P.
Dembinski15, A. Dendek34, V. Denysenko50, D. Derkach82, O. Deschamps9, F.
Desse11, F. Dettori27,f, B. Dey77, A. Di Cicco23, P. Di Nezza23, S. Didenko83,
L. Dieste Maronas46, H. Dijkstra48, V. Dobishuk52, C. Dong3, A.M. Donohoe18,
F. Dordei27, A.C. dos Reis1, L. Douglas59, A. Dovbnya51, A.G. Downes8, M.W.
Dudek35, L. Dufour48, V. Duk78, P. Durante48, J.M. Durham67, D. Dutta62, A.
Dziurda35, A. Dzyuba38, S. Easo57, U. Egede69, V. Egorychev41, S.
Eidelman43,v,†, S. Eisenhardt58, S. Ek-In49, L. Eklund86, S. Ely68, A. Ene37,
E. Epple67, S. Escher14, J. Eschle50, S. Esen50, T. Evans48, L.N. Falcao1, Y.
Fan6, B. Fang73, S. Farry60, D. Fazzini26,k, M. Féo48, A. Fernandez Prieto46,
A.D. Fernez66, F. Ferrari20,e, L. Ferreira Lopes49, F. Ferreira Rodrigues2, S.
Ferreres Sole32, M. Ferrillo50, M. Ferro-Luzzi48, S. Filippov39, R.A. Fini19,
M. Fiorini21,g, M. Firlej34, K.M. Fischer63, D.S. Fitzgerald87, C.
Fitzpatrick62, T. Fiutowski34, A. Fkiaras48, F. Fleuret12, M. Fontana13, F.
Fontanelli24,i, R. Forty48, D. Foulds-Holt55, V. Franco Lima60, M. Franco
Sevilla66, M. Frank48, E. Franzoso21, G. Frau17, C. Frei48, D.A. Friday59, J.
Fu6, Q. Fuehring15, E. Gabriel32, G. Galati19,d, A. Gallas Torreira46, D.
Galli20,e, S. Gambetta58,48, Y. Gan3, M. Gandelman2, P. Gandini25, Y. Gao5, M.
Garau27, L.M. Garcia Martin56, P. Garcia Moreno45, J. García Pardiñas26,k, B.
Garcia Plana46, F.A. Garcia Rosales12, L. Garrido45, C. Gaspar48, R.E.
Geertsema32, D. Gerick17, L.L. Gerken15, E. Gersabeck62, M. Gersabeck62, T.
Gershon56, D. Gerstel10, L. Giambastiani28, V. Gibson55, H.K. Giemza36, A.L.
Gilman63, M. Giovannetti23,q, A. Gioventù46, P. Gironella Gironell45, C.
Giugliano21,g, K. Gizdov58, E.L. Gkougkousis48, V.V. Gligorov13, C. Göbel70,
E. Golobardes85, D. Golubkov41, A. Golutvin61,83, A. Gomes1,a, S. Gomez
Fernandez45, F. Goncalves Abrantes63, M. Goncerz35, G. Gong3, P. Gorbounov41,
I.V. Gorelov40, C. Gotti26, E. Govorkova48, J.P. Grabowski17, T. Grammatico13,
L.A. Granado Cardoso48, E. Graugés45, E. Graverini49, G. Graziani22, A.
Grecu37, L.M. Greeven32, N.A. Grieser4, L. Grillo62, S. Gromov83, B.R. Gruberg
Cazon63, C. Gu3, M. Guarise21, M. Guittiere11, P. A. Günther17, E. Gushchin39,
A. Guth14, Y. Guz44, T. Gys48, T. Hadavizadeh69, G. Haefeli49, C. Haen48, J.
Haimberger48, T. Halewood-leagas60, P.M. Hamilton66, J.P. Hammerich60, Q.
Han7, X. Han17, T.H. Hancock63, E.B. Hansen62, S. Hansmann-Menzemer17, N.
Harnew63, T. Harrison60, C. Hasse48, M. Hatch48, J. He6,b, M. Hecker61, K.
Heijhoff32, K. Heinicke15, R.D.L. Henderson69,56, A.M. Hennequin48, K.
Hennessy60, L. Henry48, J. Heuel14, A. Hicheur2, D. Hill49, M. Hilton62, S.E.
Hollitt15, R. Hou7, Y. Hou8, J. Hu17, J. Hu72, W. Hu7, X. Hu3, W. Huang6, X.
Huang73, W. Hulsbergen32, R.J. Hunter56, M. Hushchyn82, D. Hutchcroft60, D.
Hynds32, P. Ibis15, M. Idzik34, D. Ilin38, P. Ilten65, A. Inglessi38, A.
Ishteev83, K. Ivshin38, R. Jacobsson48, H. Jage14, S. Jakobsen48, E. Jans32,
B.K. Jashal47, A. Jawahery66, V. Jevtic15, X. Jiang4, M. John63, D. Johnson64,
C.R. Jones55, T.P. Jones56, B. Jost48, N. Jurik48, S.H. Kalavan Kadavath34, S.
Kandybei51, Y. Kang3, M. Karacson48, M. Karpov82, J.W. Kautz65, F. Keizer48,
D.M. Keller68, M. Kenzie56, T. Ketel33, B. Khanji15, A. Kharisova84, S.
Kholodenko44, T. Kirn14, V.S. Kirsebom49, O. Kitouni64, S. Klaver32, N.
Kleijne29, K. Klimaszewski36, M.R. Kmiec36, S. Koliiev52, A. Kondybayeva83, A.
Konoplyannikov41, P. Kopciewicz34, R. Kopecna17, P. Koppenburg32, M.
Korolev40, I. Kostiuk32,52, O. Kot52, S. Kotriakhova21,38, P. Kravchenko38, L.
Kravchuk39, R.D. Krawczyk48, M. Kreps56, F. Kress61, S. Kretzschmar14, P.
Krokovny43,v, W. Krupa34, W. Krzemien36, J. Kubat17, M. Kucharczyk35, V.
Kudryavtsev43,v, H.S. Kuindersma32,33, G.J. Kunde67, T. Kvaratskheliya41, D.
Lacarrere48, G. Lafferty62, A. Lai27, A. Lampis27, D. Lancierini50, J.J.
Lane62, R. Lane54, G. Lanfranchi23, C. Langenbruch14, J. Langer15, O.
Lantwin83, T. Latham56, F. Lazzari29,r, R. Le Gac10, S.H. Lee87, R. Lefèvre9,
A. Leflat40, S. Legotin83, O. Leroy10, T. Lesiak35, B. Leverington17, H. Li72,
P. Li17, S. Li7, Y. Li4, Y. Li4, Z. Li68, X. Liang68, T. Lin61, R. Lindner48,
V. Lisovskyi15, R. Litvinov27, G. Liu72, H. Liu6, Q. Liu6, S. Liu4, A. Lobo
Salvia45, A. Loi27, J. Lomba Castro46, I. Longstaff59, J.H. Lopes2, S. López
Soliño46, G.H. Lovell55, Y. Lu4, C. Lucarelli22,h, D. Lucchesi28,m, S.
Luchuk39, M. Lucio Martinez32, V. Lukashenko32,52, Y. Luo3, A. Lupato62, E.
Luppi21,g, O. Lupton56, A. Lusiani29,n, X. Lyu6, L. Ma4, R. Ma6, S.
Maccolini20,e, F. Machefert11, F. Maciuc37, V. Macko49, P. Mackowiak15, S.
Maddrell-Mander54, O. Madejczyk34, L.R. Madhan Mohan54, O. Maev38, A.
Maevskiy82, M.W. Majewski34, J.J. Malczewski35, S. Malde63, B. Malecki48, A.
Malinin81, T. Maltsev43,v, H. Malygina17, G. Manca27,f, G. Mancinelli10, D.
Manuzzi20,e, D. Marangotto25,j, J. Maratas9,t, J.F. Marchand8, U. Marconi20,
S. Mariani22,h, C. Marin Benito48, M. Marinangeli49, J. Marks17, A.M.
Marshall54, P.J. Marshall60, G. Martelli78, G. Martellotti30, L.
Martinazzoli48,k, M. Martinelli26,k, D. Martinez Santos46, F. Martinez
Vidal47, A. Massafferri1, M. Materok14, R. Matev48, A. Mathad50, V.
Matiunin41, C. Matteuzzi26, K.R. Mattioli87, A. Mauri32, E. Maurice12, J.
Mauricio45, M. Mazurek48, M. McCann61, L. Mcconnell18, T.H. Mcgrath62, N.T.
Mchugh59, A. McNab62, R. McNulty18, J.V. Mead60, B. Meadows65, G. Meier15, D.
Melnychuk36, S. Meloni26,k, M. Merk32,80, A. Merli25,j, L. Meyer Garcia2, M.
Mikhasenko75,c, D.A. Milanes74, E. Millard56, M. Milovanovic48, M.-N. Minard8,
A. Minotti26,k, L. Minzoni21,g, S.E. Mitchell58, B. Mitreska62, D.S. Mitzel15,
A. Mödden 15, R.A. Mohammed63, R.D. Moise61, S. Mokhnenko82, T. Mombächer46,
I.A. Monroy74, S. Monteil9, M. Morandin28, G. Morello23, M.J. Morello29,n, J.
Moron34, A.B. Morris75, A.G. Morris56, R. Mountain68, H. Mu3, F. Muheim58,48,
M. Mulder79, D. Müller48, K. Müller50, C.H. Murphy63, D. Murray62, R. Murta61,
P. Muzzetto27, P. Naik54, T. Nakada49, R. Nandakumar57, T. Nanut48, I.
Nasteva2, M. Needham58, N. Neri25,j, S. Neubert75, N. Neufeld48, R.
Newcombe61, E.M. Niel11, S. Nieswand14, N. Nikitin40, N.S. Nolte64, C.
Normand8, C. Nunez87, A. Oblakowska-Mucha34, V. Obraztsov44, T. Oeser14, D.P.
O’Hanlon54, S. Okamura21, R. Oldeman27,f, F. Oliva58, M.E. Olivares68, C.J.G.
Onderwater79, R.H. O’Neil58, J.M. Otalora Goicochea2, T. Ovsiannikova41, P.
Owen50, A. Oyanguren47, K.O. Padeken75, B. Pagare56, P.R. Pais48, T. Pajero63,
A. Palano19, M. Palutan23, Y. Pan62, G. Panshin84, A. Papanestis57, M.
Pappagallo19,d, L.L. Pappalardo21,g, C. Pappenheimer65, W. Parker66, C.
Parkes62, B. Passalacqua21, G. Passaleva22, A. Pastore19, M. Patel61, C.
Patrignani20,e, C.J. Pawley80, A. Pearce48,57, A. Pellegrino32, M. Pepe
Altarelli48, S. Perazzini20, D. Pereima41, A. Pereiro Castro46, P. Perret9, M.
Petric59,48, K. Petridis54, A. Petrolini24,i, A. Petrov81, S. Petrucci58, M.
Petruzzo25, T.T.H. Pham68, A. Philippov42, R. Piandani6, L. Pica29,n, M.
Piccini78, B. Pietrzyk8, G. Pietrzyk49, M. Pili63, D. Pinci30, F. Pisani48, M.
Pizzichemi26,48,k, Resmi P.K10, V. Placinta37, J. Plews53, M. Plo Casasus46,
F. Polci13, M. Poli Lener23, M. Poliakova68, A. Poluektov10, N. Polukhina83,u,
I. Polyakov68, E. Polycarpo2, S. Ponce48, D. Popov6,48, S. Popov42, S.
Poslavskii44, K. Prasanth35, L. Promberger48, C. Prouve46, V. Pugatch52, V.
Puill11, G. Punzi29,o, H. Qi3, W. Qian6, N. Qin3, R. Quagliani49, N.V. Raab18,
R.I. Rabadan Trejo6, B. Rachwal34, J.H. Rademacker54, M. Rama29, M. Ramos
Pernas56, M.S. Rangel2, F. Ratnikov42,82, G. Raven33, M. Reboud8, F. Redi49,
F. Reiss62, C. Remon Alepuz47, Z. Ren3, V. Renaudin63, R. Ribatti29, A.M.
Ricci27, S. Ricciardi57, K. Rinnert60, P. Robbe11, G. Robertson58, A.B.
Rodrigues49, E. Rodrigues60, J.A. Rodriguez Lopez74, E.R.R. Rodriguez
Rodriguez46, A. Rollings63, P. Roloff48, V. Romanovskiy44, M. Romero Lamas46,
A. Romero Vidal46, J.D. Roth87, M. Rotondo23, M.S. Rudolph68, T. Ruf48, R.A.
Ruiz Fernandez46, J. Ruiz Vidal47, A. Ryzhikov82, J. Ryzka34, J.J. Saborido
Silva46, N. Sagidova38, N. Sahoo56, B. Saitta27,f, M. Salomoni48, C. Sanchez
Gras32, R. Santacesaria30, C. Santamarina Rios46, M. Santimaria23, E.
Santovetti31,q, D. Saranin83, G. Sarpis14, M. Sarpis75, A. Sarti30, C.
Satriano30,p, A. Satta31, M. Saur15, D. Savrina41,40, H. Sazak9, L.G.
Scantlebury Smead63, A. Scarabotto13, S. Schael14, S. Scherl60, M. Schiller59,
H. Schindler48, M. Schmelling16, B. Schmidt48, S. Schmitt14, O. Schneider49,
A. Schopper48, M. Schubiger32, S. Schulte49, M.H. Schune11, R. Schwemmer48, B.
Sciascia23,48, S. Sellam46, A. Semennikov41, M. Senghi Soares33, A. Sergi24,i,
N. Serra50, L. Sestini28, A. Seuthe15, Y. Shang5, D.M. Shangase87, M.
Shapkin44, I. Shchemerov83, L. Shchutska49, T. Shears60, L. Shekhtman43,v, Z.
Shen5, S. Sheng4, V. Shevchenko81, E.B. Shields26,k, Y. Shimizu11, E.
Shmanin83, J.D. Shupperd68, B.G. Siddi21, R. Silva Coutinho50, G. Simi28, S.
Simone19,d, N. Skidmore62, T. Skwarnicki68, M.W. Slater53, I. Slazyk21,g, J.C.
Smallwood63, J.G. Smeaton55, A. Smetkina41, E. Smith50, M. Smith61, A.
Snoch32, L. Soares Lavra9, M.D. Sokoloff65, F.J.P. Soler59, A. Solovev38, I.
Solovyev38, F.L. Souza De Almeida2, B. Souza De Paula2, B. Spaan15, E. Spadaro
Norella25,j, P. Spradlin59, F. Stagni48, M. Stahl65, S. Stahl48, S.
Stanislaus63, O. Steinkamp50,83, O. Stenyakin44, H. Stevens15, S.
Stone68,48,†, D. Strekalina83, F. Suljik63, J. Sun27, L. Sun73, Y. Sun66, P.
Svihra62, P.N. Swallow53, K. Swientek34, A. Szabelski36, T. Szumlak34, M.
Szymanski48, S. Taneja62, A.R. Tanner54, M.D. Tat63, A. Terentev83, F.
Teubert48, E. Thomas48, D.J.D. Thompson53, K.A. Thomson60, H. Tilquin61, V.
Tisserand9, S. T’Jampens8, M. Tobin4, L. Tomassetti21,g, X. Tong5, D. Torres
Machado1, D.Y. Tou13, E. Trifonova83, S.M. Trilov54, C. Trippl49, G. Tuci6, A.
Tully49, N. Tuning32,48, A. Ukleja36,48, D.J. Unverzagt17, E. Ursov83, A.
Usachov32, A. Ustyuzhanin42,82, U. Uwer17, A. Vagner84, V. Vagnoni20, A.
Valassi48, G. Valenti20, N. Valls Canudas85, M. van Beuzekom32, M. Van Dijk49,
H. Van Hecke67, E. van Herwijnen83, M. van Veghel79, R. Vazquez Gomez45, P.
Vazquez Regueiro46, C. Vázquez Sierra48, S. Vecchi21, J.J. Velthuis54, M.
Veltri22,s, A. Venkateswaran68, M. Veronesi32, M. Vesterinen56, D. Vieira65,
M. Vieites Diaz49, H. Viemann76, X. Vilasis-Cardona85, E. Vilella Figueras60,
A. Villa20, P. Vincent13, F.C. Volle11, D. Vom Bruch10, A. Vorobyev38,†, V.
Vorobyev43,v, N. Voropaev38, K. Vos80, R. Waldi17, J. Walsh29, C. Wang17, J.
Wang5, J. Wang4, J. Wang3, J. Wang73, M. Wang3, R. Wang54, Y. Wang7, Z.
Wang50, Z. Wang3, Z. Wang6, J.A. Ward56,69, N.K. Watson53, S.G. Weber13, D.
Websdale61, C. Weisser64, B.D.C. Westhenry54, D.J. White62, M. Whitehead54,
A.R. Wiederhold56, D. Wiedner15, G. Wilkinson63, M. Wilkinson68, I.
Williams55, M. Williams64, M.R.J. Williams58, F.F. Wilson57, W. Wislicki36, M.
Witek35, L. Witola17, G. Wormser11, S.A. Wotton55, H. Wu68, K. Wyllie48, Z.
Xiang6, D. Xiao7, Y. Xie7, A. Xu5, J. Xu6, L. Xu3, M. Xu56, Q. Xu6, Z. Xu9, Z.
Xu6, D. Yang3, S. Yang6, Y. Yang6, Z. Yang5, Z. Yang66, Y. Yao68, L.E.
Yeomans60, H. Yin7, J. Yu71, X. Yuan68, O. Yushchenko44, E. Zaffaroni49, M.
Zavertyaev16,u, M. Zdybal35, O. Zenaiev48, M. Zeng3, D. Zhang7, L. Zhang3, S.
Zhang71, S. Zhang5, Y. Zhang5, Y. Zhang63, A. Zharkova83, A. Zhelezov17, Y.
Zheng6, T. Zhou5, X. Zhou6, Y. Zhou6, V. Zhovkovska11, X. Zhu3, X. Zhu7, Z.
Zhu6, V. Zhukov14,40, J.B. Zonneveld58, Q. Zou4, S. Zucchelli20,e, D.
Zuliani28, G. Zunica62.

1Centro Brasileiro de Pesquisas Físicas (CBPF), Rio de Janeiro, Brazil
2Universidade Federal do Rio de Janeiro (UFRJ), Rio de Janeiro, Brazil
3Center for High Energy Physics, Tsinghua University, Beijing, China
4Institute Of High Energy Physics (IHEP), Beijing, China
5School of Physics State Key Laboratory of Nuclear Physics and Technology,
Peking University, Beijing, China
6University of Chinese Academy of Sciences, Beijing, China
7Institute of Particle Physics, Central China Normal University, Wuhan, Hubei,
China
8Univ. Savoie Mont Blanc, CNRS, IN2P3-LAPP, Annecy, France
9Université Clermont Auvergne, CNRS/IN2P3, LPC, Clermont-Ferrand, France
10Aix Marseille Univ, CNRS/IN2P3, CPPM, Marseille, France
11Université Paris-Saclay, CNRS/IN2P3, IJCLab, Orsay, France
12Laboratoire Leprince-Ringuet, CNRS/IN2P3, Ecole Polytechnique, Institut
Polytechnique de Paris, Palaiseau, France
13LPNHE, Sorbonne Université, Paris Diderot Sorbonne Paris Cité, CNRS/IN2P3,
Paris, France
14I. Physikalisches Institut, RWTH Aachen University, Aachen, Germany
15Fakultät Physik, Technische Universität Dortmund, Dortmund, Germany
16Max-Planck-Institut für Kernphysik (MPIK), Heidelberg, Germany
17Physikalisches Institut, Ruprecht-Karls-Universität Heidelberg, Heidelberg,
Germany
18School of Physics, University College Dublin, Dublin, Ireland
19INFN Sezione di Bari, Bari, Italy
20INFN Sezione di Bologna, Bologna, Italy
21INFN Sezione di Ferrara, Ferrara, Italy
22INFN Sezione di Firenze, Firenze, Italy
23INFN Laboratori Nazionali di Frascati, Frascati, Italy
24INFN Sezione di Genova, Genova, Italy
25INFN Sezione di Milano, Milano, Italy
26INFN Sezione di Milano-Bicocca, Milano, Italy
27INFN Sezione di Cagliari, Monserrato, Italy
28Universita degli Studi di Padova, Universita e INFN, Padova, Padova, Italy
29INFN Sezione di Pisa, Pisa, Italy
30INFN Sezione di Roma La Sapienza, Roma, Italy
31INFN Sezione di Roma Tor Vergata, Roma, Italy
32Nikhef National Institute for Subatomic Physics, Amsterdam, Netherlands
33Nikhef National Institute for Subatomic Physics and VU University Amsterdam,
Amsterdam, Netherlands
34AGH - University of Science and Technology, Faculty of Physics and Applied
Computer Science, Kraków, Poland
35Henryk Niewodniczanski Institute of Nuclear Physics Polish Academy of
Sciences, Kraków, Poland
36National Center for Nuclear Research (NCBJ), Warsaw, Poland
37Horia Hulubei National Institute of Physics and Nuclear Engineering,
Bucharest-Magurele, Romania
38Petersburg Nuclear Physics Institute NRC Kurchatov Institute (PNPI NRC KI),
Gatchina, Russia
39Institute for Nuclear Research of the Russian Academy of Sciences (INR RAS),
Moscow, Russia
40Institute of Nuclear Physics, Moscow State University (SINP MSU), Moscow,
Russia
41Institute of Theoretical and Experimental Physics NRC Kurchatov Institute
(ITEP NRC KI), Moscow, Russia
42Yandex School of Data Analysis, Moscow, Russia
43Budker Institute of Nuclear Physics (SB RAS), Novosibirsk, Russia
44Institute for High Energy Physics NRC Kurchatov Institute (IHEP NRC KI),
Protvino, Russia, Protvino, Russia
45ICCUB, Universitat de Barcelona, Barcelona, Spain
46Instituto Galego de Física de Altas Enerxías (IGFAE), Universidade de
Santiago de Compostela, Santiago de Compostela, Spain
47Instituto de Fisica Corpuscular, Centro Mixto Universidad de Valencia -
CSIC, Valencia, Spain
48European Organization for Nuclear Research (CERN), Geneva, Switzerland
49Institute of Physics, Ecole Polytechnique Fédérale de Lausanne (EPFL),
Lausanne, Switzerland
50Physik-Institut, Universität Zürich, Zürich, Switzerland
51NSC Kharkiv Institute of Physics and Technology (NSC KIPT), Kharkiv, Ukraine
52Institute for Nuclear Research of the National Academy of Sciences (KINR),
Kyiv, Ukraine
53University of Birmingham, Birmingham, United Kingdom
54H.H. Wills Physics Laboratory, University of Bristol, Bristol, United
Kingdom
55Cavendish Laboratory, University of Cambridge, Cambridge, United Kingdom
56Department of Physics, University of Warwick, Coventry, United Kingdom
57STFC Rutherford Appleton Laboratory, Didcot, United Kingdom
58School of Physics and Astronomy, University of Edinburgh, Edinburgh, United
Kingdom
59School of Physics and Astronomy, University of Glasgow, Glasgow, United
Kingdom
60Oliver Lodge Laboratory, University of Liverpool, Liverpool, United Kingdom
61Imperial College London, London, United Kingdom
62Department of Physics and Astronomy, University of Manchester, Manchester,
United Kingdom
63Department of Physics, University of Oxford, Oxford, United Kingdom
64Massachusetts Institute of Technology, Cambridge, MA, United States
65University of Cincinnati, Cincinnati, OH, United States
66University of Maryland, College Park, MD, United States
67Los Alamos National Laboratory (LANL), Los Alamos, United States
68Syracuse University, Syracuse, NY, United States
69School of Physics and Astronomy, Monash University, Melbourne, Australia,
associated to 56
70Pontifícia Universidade Católica do Rio de Janeiro (PUC-Rio), Rio de
Janeiro, Brazil, associated to 2
71Physics and Micro Electronic College, Hunan University, Changsha City,
China, associated to 7
72Guangdong Provincial Key Laboratory of Nuclear Science, Guangdong-Hong Kong
Joint Laboratory of Quantum Matter, Institute of Quantum Matter, South China
Normal University, Guangzhou, China, associated to 3
73School of Physics and Technology, Wuhan University, Wuhan, China, associated
to 3
74Departamento de Fisica , Universidad Nacional de Colombia, Bogota, Colombia,
associated to 13
75Universität Bonn - Helmholtz-Institut für Strahlen und Kernphysik, Bonn,
Germany, associated to 17
76Institut für Physik, Universität Rostock, Rostock, Germany, associated to 17
77Eotvos Lorand University, Budapest, Hungary, associated to 48
78INFN Sezione di Perugia, Perugia, Italy, associated to 21
79Van Swinderen Institute, University of Groningen, Groningen, Netherlands,
associated to 32
80Universiteit Maastricht, Maastricht, Netherlands, associated to 32
81National Research Centre Kurchatov Institute, Moscow, Russia, associated to
41
82National Research University Higher School of Economics, Moscow, Russia,
associated to 42
83National University of Science and Technology “MISIS”, Moscow, Russia,
associated to 41
84National Research Tomsk Polytechnic University, Tomsk, Russia, associated to
41
85DS4DS, La Salle, Universitat Ramon Llull, Barcelona, Spain, associated to 45
86Department of Physics and Astronomy, Uppsala University, Uppsala, Sweden,
associated to 59
87University of Michigan, Ann Arbor, United States, associated to 68

aUniversidade Federal do Triângulo Mineiro (UFTM), Uberaba-MG, Brazil
bHangzhou Institute for Advanced Study, UCAS, Hangzhou, China
cExcellence Cluster ORIGINS, Munich, Germany
dUniversità di Bari, Bari, Italy
eUniversità di Bologna, Bologna, Italy
fUniversità di Cagliari, Cagliari, Italy
gUniversità di Ferrara, Ferrara, Italy
hUniversità di Firenze, Firenze, Italy
iUniversità di Genova, Genova, Italy
jUniversità degli Studi di Milano, Milano, Italy
kUniversità di Milano Bicocca, Milano, Italy
lUniversità di Modena e Reggio Emilia, Modena, Italy
mUniversità di Padova, Padova, Italy
nScuola Normale Superiore, Pisa, Italy
oUniversità di Pisa, Pisa, Italy
pUniversità della Basilicata, Potenza, Italy
qUniversità di Roma Tor Vergata, Roma, Italy
rUniversità di Siena, Siena, Italy
sUniversità di Urbino, Urbino, Italy
tMSU - Iligan Institute of Technology (MSU-IIT), Iligan, Philippines
uP.N. Lebedev Physical Institute, Russian Academy of Science (LPI RAS),
Moscow, Russia
vNovosibirsk State University, Novosibirsk, Russia
†Deceased


Same as section 5.3.4.3 but with $i$ and $k$ exchanged, and with $j$ and
$\ell$ exchanged.

#####  5.3.4.5 $\ell=j$ and $i$ and $k$ are distinct

The clusters are $ij$ and $kj$. So, by the same analysis from section 5.3.4.1,
we get

$\displaystyle\sum_{n=0}^{\infty}E[X_{i,j}^{(n)}X_{k,j}^{(n)}]z^{n}$
$\displaystyle=\frac{1}{1-z}-F(z,0,1)-F(z,1,0)+F(z,0,0)$
$\displaystyle=\frac{1}{1-z}-\frac{1}{1-z+P_{i}P_{j}z^{2}}-\frac{1}{1-z+P_{k}P_{j}z^{2}}$
$\displaystyle\qquad{}+\frac{1}{1-z+P_{i}P_{j}z^{2}+P_{k}P_{j}z^{2}}$

Exactly as in section 5.3.4.1 above:

The partial fraction decomposition for the second term is given in Table 1A.

The partial fraction decomposition for the third term is given in Table 1A,
using $k$ instead of $i$.

The partial fraction decomposition for the fourth term is given in Table 1C,
using $j$ instead of $\ell$.

#####  5.3.4.6 $i=k$ and $j=\ell$ are distinct

In this case we have $X_{i,j}^{(n)}X_{k,\ell}^{(n)}=X_{i,j}^{(n)}$, so we get
$E[X_{i,j}^{(n)}X_{k,\ell}^{(n)}]=E[X_{i,j}^{(n)}]$, which we already handled
in Lemma 5.2.

#####  5.3.4.7 $i=\ell$ and $j=k$ are distinct

The clusters each have the form $ijiji\ldots$ or $jijij\ldots$. So
$\xi(z,t,u)$ of the set of clusters
$C=\\{ij,iji,ijij,ijiji\ldots,ji,jij,jiji,jijij\ldots\\}$ becomes

$\xi(z,t,u)=\frac{P_{i}P_{j}tz^{2}}{1-P_{i}P_{j}tuz^{2}}+\frac{P_{i}^{2}P_{j}tuz^{3}}{1-P_{i}P_{j}tuz^{2}}+\frac{P_{j}P_{i}uz^{2}}{1-P_{i}P_{j}tuz^{2}}+\frac{P_{i}P_{j}^{2}tuz^{3}}{1-P_{i}P_{j}tuz^{2}},$

and $F(z,x,y)=1/(1-z-\xi(z,x-1,y-1))$. It follows that the probability
generating function of the words with _at least one occurrence of_ $ij$ and
_at least one occurrence of_ $k\ell$ is

$\displaystyle\sum_{n=0}^{\infty}E[X_{i,j}^{(n)}X_{k,\ell}^{(n)}]z^{n}$
$\displaystyle=\frac{1}{1-z}-F(z,0,1)-F(z,1,0)+F(z,0,0)$
$\displaystyle=\frac{1}{1-z}-\frac{1}{1-z+P_{i}P_{j}z^{2}}-\frac{1}{1-z+P_{j}P_{i}z^{2}}$
$\displaystyle\qquad{}+\frac{1}{1-z+\frac{P_{i}P_{j}z^{2}}{1-P_{i}P_{j}z^{2}}-\frac{P_{i}^{2}P_{j}z^{3}}{1-P_{i}P_{j}z^{2}}+\frac{P_{j}P_{i}z^{2}}{1-P_{i}P_{j}z^{2}}-\frac{P_{i}P_{j}^{2}z^{3}}{1-P_{i}P_{j}z^{2}}}$

The partial fraction decomposition for the second and for the third term is
given in Table 1A.

The partial fraction decomposition for the fourth term is given in Table 2H.

As mentioned immediately after Theorem 5.4, we do not wrap all of the analysis
from Section 5.3 into a theorem (because it would be very lengthy), but we
have precisely analyzed every aspect that is needed for exactly characterizing
the second moment $E[(X^{(n)})^{2}]$.

###### Acknowledgements.

We would like to thank B. Pittel for providing the Poisson distribution of
$X_{i,j}(m)$, and B. Salvy for suggesting the use of the Maple package gfun.
Furthermore we would like to thank an anonymous referee, whose suggestions led
to substantial improvements of the paper. In particular we want to thank for
pointing out to us a connection to Stirling numbers, that is stated in Remark
3.3. M.D. Ward’s research is supported by National Science Foundation (NSF)
grants 0939370, 1246818, 2005632, 2123321, 2118329, and 2235473, by the
Foundation for Food and Agriculture Research (FFAR) grant 534662, by the
National Institute of Food and Agriculture (NIFA) grants 2019-67032-29077,
2020-70003-32299, 2021-38420-34943, and 2022-67021-37022, by the Society Of
Actuaries grant 19111857, by Cummins Inc., by Gro Master, by Lilly Endowment,
and by Sandia National Laboratories.

## References

  * Archibald et al. (2021) M. Archibald, A. Blecher, C. Brennan, A. Knopfmacher, S. Wagner, and M. Ward.  The number of distinct adjacent pairs in geometrically distributed words.  _Discrete Mathematics and Theoretical Computer Science_ , 22(4), 2021.
  * Bassino et al. (2012) F. Bassino, J. Clément, and P. Nicodème.  Counting occurrences for a finite set of words: Combinatorial methods.  _ACM Transactions on Algorithms_ , 8(3), 2012.  Article 31.
  * Flajolet and Sedgewick (2009) P. Flajolet and R. Sedgewick.  _Analytic Combinatorics_.  Cambridge, 2009.
  * Flajolet et al. (1995) P. Flajolet, X. Gourdon, and P. Dumas.  Mellin transforms and asymptotics: Harmonic sums.  _Theoretical Computer Science_ , 144:3–58, 1995.
  * Hitczenko and Louchard (2001) P. Hitczenko and G. Louchard.  Distinctness of compositions of an integer: A probabilistic analysis.  _Random Structures and Algorithms_ , 19(3–4):407–437, 2001.
  * Louchard and Prodinger (2006) G. Louchard and H. Prodinger.  Asymptotics of the moments of extreme-value related distribution functions.  _Algorithmica_ , 46:431–467, 2006.
  * Louchard et al. (2005) G. Louchard, H. Prodinger, and M. Ward.  The number of distinct values of some multiplicity in sequences of geometrically distributed random variables.  _Discrete Mathematics and Theoretical Computer Science_ , AD:231–256, 2005.  Proceedings of the 2005 International Conference on Analysis of Algorithms.
  * Louchard et al. (2023) G. Louchard, W. Schachinger, and M. Ward.  The number of distinct values of some multiplicity in sequences of geometrically distributed random variables.  _arXiv_ , 2023.  Extended preprint https://arxiv.org/abs/2203.14773.
  * McDiarmid (1992) C. McDiarmid.  On a correlation inequality of Farr.  _Combinatorics, Probability and Computing_ , 1:157–160, 1992.
  * (10) B. Pittel.  Technical report. private communication.
  * Rudin (1976) W. Rudin.  _Principles of Mathematical Analysis, 3rd ed._ McGraw-Hill, 1976.
  * (12) B. Salvy.  Private communication.
  * Salvy and Zimmermann (1994) B. Salvy and P. Zimmermann.  Gfun: A Maple package for the manipulation of generating and holonomic functions in one variable.  _ACM Transactions on Mathematical Software_ , 20(2):163–177, 1994.
  * Seneta (1981) E. Seneta.  _Non-negative Matrices and Markov Chains, 2nd ed._ Springer, 1981.
  * (15) N. J. A. Sloane.  The on-line encyclopedia of integer sequences.  http://oeis.org. Sequence A028246.

##  Appendix A Some Mellin transforms

To keep the paper self contained we give here a short outline on how to use
Mellin transforms to obtain asymptotic expansions. The reader seeking more
detail is referred to Flajolet et al. (1995) for a nice exposition.
Subsections A.1, A.2, and A.3 are devoted to asymptotic equivalents of three
sums that play a crucial role in our paper.

The Mellin transform $f^{*}(s)$ of $f(x)$, also denoted
$\mathcal{M}\left[f(x);s\right]$, is given by

$f^{*}(s)=\int_{0}^{\infty}f(x)x^{s-1}dx.$

The interior of the set of $s$ for which the integral converges is an open
strip $\langle a,b\rangle:=\\{s\in\mathbb{C}:a<\Re s<b\\}$, called the
_fundamental strip_ , with $a,b$ depending on how $f$ behaves at $0$ and
$\infty$. For example, we have $\mathcal{M}\left[e^{-x};s\right]=\Gamma(s)$,
with fundamental strip $\langle 0,\infty\rangle$, and
$\mathcal{M}\left[1-e^{-x};s\right]=-\Gamma(s)$, with fundamental strip
$\langle-1,0\rangle$. When computing the Mellin transform of so called
_harmonic sums_ , the rescaling rule turns out to be very useful:

$\mathcal{M}\left[\sum_{k}\lambda_{k}f(\mu_{k}x);s\right]=\sum_{k}\frac{\lambda_{k}}{\mu_{k}^{s}}\cdot
f^{*}(s).$

In the case that $f^{*}(s)$ can be meromorphically continued to a strip
$\langle a,\bar{b}\rangle$ with $\bar{b}>b$, information on the poles of
$f^{*}(s)$ leads to asymptotic properties of $f(n),n\rightarrow\infty$. This
is called the fundamental correspondence. In particular, if there is a pole

$\frac{1}{(s-\xi)^{k+1}}$

of $f^{*}(s)$ at $\xi=\sigma+\mathbf{i}t$ to the right of the fundamental
strip, then this pole will contribute the term

$-\frac{(-1)^{k}}{k!}\ln(n)^{k}n^{-\sigma}e^{-\mathbf{i}t\ln(n)},$

which is precisely the residue of $\frac{-n^{-s}}{(s-\xi)^{k+1}}$ at $s=\xi$,
to an asymptotic expansion of $f(n)$ at $\infty$. Justification comes from
residue calculus: If $f$ is smooth enough, the inverse transform applies to
yield
$f(n)=\frac{1}{2\pi\mathbf{i}}\int_{c-\mathbf{i}\infty}^{c+\mathbf{i}\infty}n^{-s}f^{*}(s)ds$
with $a<c<b$. If for $\beta<\bar{b}$ the set of poles $z$ in $\langle
a,\bar{b}\rangle$ satisfying $\Re z<\beta$ is denoted $\mathcal{S}_{\beta}$,
and there is no pole with real part $\beta$, we have

$f(n)=\frac{1}{2\pi\mathbf{i}}\int_{\beta-\mathbf{i}\infty}^{\beta+\mathbf{i}\infty}n^{-s}f^{*}(s)ds-\sum_{z\in\mathcal{S}_{\beta}}\textup{Res}(f^{*}(s)n^{-s})_{s=z}=-\sum_{z\in\mathcal{S}_{\beta}}\textup{Res}(f^{*}(s)n^{-s})_{s=z}+\mathcal{O}(n^{-\beta}),$

with the integral being $\mathcal{O}(n^{-\beta})$, provided that $f^{*}(s)$
decreases fast enough for $s=\beta+\mathbf{i}t$ and $|t|\to\infty$. Equality
of left and right hand side is established by using a sequence of contours
$\rho_{k}$ being the boundaries of rectangles $\\{z\in\mathbb{C}:c\leq\Re
z\leq\beta,-h_{k}\leq\Im z\leq h_{k}\\}$ with $h_{k}\to\infty$, verifying that
$f^{*}(s)$ decreases fast enough on the horizontal segments of $\rho_{k}$, as
$k\to\infty$, and applying residue calculus.

Here is an illustration of fast enough decrease. $\Gamma(s)$ decreases
exponentially in the direction $\mathbf{i}\infty$:

$|\Gamma(\sigma+\mathbf{i}t)|\sim\sqrt{2\pi}|t|^{\sigma-1/2}e^{-\pi|t|/2}.$

Also, similarly fast decrease can be observed for all other transforms we
encounter.

In the following, recall the notations

$L:=\ln\tfrac{1}{q}\mbox{ and }\chi:=\tfrac{2\mathbf{i}\pi}{L}.$

### A.1

Let

$G(n):=\sum_{i\geq 0}\left(1-e^{-nq^{2i}}\right).$

The Mellin transform of this sum is $G^{*}(s)=-\dfrac{1}{1-q^{-2s}}\Gamma(s)$,
with fundamental strip $\langle-1,0\rangle$, to the right of which the
meromorphic extension of $G^{*}(s)$ has poles at $s=0$ and
$s=\frac{\ell\chi}{2}$ for $\ell\in\mbox{$\mathbb{Z}$}\setminus\\{0\\}$, with
singular expansions

$G^{*}(s)\asymp\frac{1}{2Ls^{2}}-\left[\frac{\gamma}{2L}+\frac{1}{2}\right]\frac{1}{s},\quad\mbox{
and
}G^{*}(s)\asymp\frac{1}{2L}\frac{\Gamma\big{(}\frac{\ell\chi}{2}\big{)}}{\big{(}s-\frac{\ell\chi}{2}\big{)}},\mbox{
for }\ell\in\mbox{$\mathbb{Z}$}\setminus\\{0\\}.$

Noting that there are no other singularities to the right of
$\langle-1,0\rangle$, the error term in the following expansion can be chosen
$\mathcal{O}(n^{-\beta})$ with any fixed $\beta>0$.

$G(n)\sim\frac{1}{2L}\ln(n)+\left[\frac{\gamma}{2L}+\frac{1}{2}\right]-\frac{1}{2L}\sum_{\ell\in\mbox{$\mathbb{Z}$}\setminus\\{0\\}}\Gamma\Big{(}\frac{\ell\chi}{2}\Big{)}n^{-\ell\chi/2}.$

### A.2

Let

$\tilde{G}(n):=\sum_{i,j\geq 0}\left(1-e^{-nq^{i+j}}\right)=\sum_{k\geq
0}(k+1)\left(1-e^{-nq^{k}}\right).$

Here we have $\tilde{G}^{*}(s)=-\dfrac{1}{(1-q^{-s})^{2}}\Gamma(s)$, with
fundamental strip $\langle-1,0\rangle$, and poles at $s=0$ and $s=\ell\chi$
for $\ell\in\mbox{$\mathbb{Z}$}\setminus\\{0\\}$, with singular expansions

$\tilde{G}^{*}(s)\asymp-\frac{1}{L^{2}s^{3}}+\left[\frac{\gamma}{L^{2}}+\frac{1}{L}\right]\frac{1}{s^{2}}-\left[\frac{\pi^{2}+6\gamma^{2}}{12L^{2}}+\frac{5}{12}+\frac{\gamma}{L}\right]\frac{1}{s},$

and

$\tilde{G}^{*}(s)\asymp-\frac{\Gamma(\ell\chi)}{L^{2}(s-\ell\chi)^{2}}-\frac{\Gamma^{\prime}(\ell\chi)-L\Gamma(\ell\chi)}{L^{2}(s-\ell\chi)},\mbox{
for }\ell\in\mbox{$\mathbb{Z}$}\setminus\\{0\\},$

leading to

$\tilde{G}(n)\sim\frac{\ln(n)^{2}}{2L^{2}}+\left[\frac{\gamma}{L^{2}}+\frac{1}{L}\right]\ln(n)+\left[\frac{\pi^{2}+6\gamma^{2}}{12L^{2}}+\frac{5}{12}+\frac{\gamma}{L}\right]+\frac{1}{L^{2}}\\!\\!\\!\sum_{\ell\in\mbox{$\mathbb{Z}$}\setminus\\{0\\}}\\!\\!\\!\\!\\!\\!\left[\Gamma^{\prime}\left(\ell\chi\right)-(\ln(n)+L)\Gamma\left(\ell\chi\right)\right]n^{-\ell\chi},$

again with error term $\mathcal{O}(n^{-\beta})$ with any fixed $\beta>0$.

### A.3

Set

$\hat{G}(n):=\sum_{i,j,k}(e^{nP_{i}P_{j}P_{k}}-1)e^{-nP_{i}P_{j}-nP_{j}P_{k}}.$

This leads to the Mellin transform, with fundamental strip
$\langle-1,0\rangle$,

$\displaystyle\hat{G}^{*}(s)=$
$\displaystyle\sum_{i,j,k}\int_{0}^{\infty}(e^{xP_{i}P_{j}P_{k}}-1)e^{-xP_{i}P_{j}-xP_{j}P_{k}}x^{s-1}dx$
$\displaystyle=$
$\displaystyle\sum_{j}P_{j}^{-s}\sum_{i,k}\int_{0}^{\infty}\left[e^{-y(P_{i}+P_{k}-P_{i}P_{k})}-e^{-y(P_{i}+P_{k})}\right]y^{s-1}dy$
$\displaystyle=$
$\displaystyle\frac{q^{s}}{p^{s}(q^{s}-1)}\Gamma(s)\sum_{i,k}\left[(P_{i}+P_{k}-P_{i}P_{k})^{-s}-(P_{i}+P_{k})^{-s}\right]$
$\displaystyle=$
$\displaystyle\left(\frac{q}{p}\right)^{2s}\frac{\Gamma(s)}{q^{s}-1}F_{1}(s),$

where

$\displaystyle F_{1}(s)$
$\displaystyle=\sum_{i,k}\left[(q^{i}+q^{k}-pq^{i+k-1})^{-s}-(q^{i}+q^{k})^{-s}\right]$
$\displaystyle=\sum_{i\geq 1}q^{-is}\left[(2-P_{i})^{-s}-2^{-s}+2\sum_{j\geq
1}\left[(1+q^{j}-pq^{i+j-1})^{-s}-(1+q^{j})^{-s}\right]\right].$

Note that $F_{1}(s)$, being a general Dirichlet series in the variable $-s$,
is analytic at least for $\sigma=\Re s<1$, since, using the Mean Value
Theorem, we have

$\left|(1+q^{j}-pq^{i+j-1})^{-\sigma}-(1+q^{j})^{-\sigma}\right|\leq|\sigma|\frac{pq^{i+j-1}}{(1+q^{j})^{1+\sigma}},$

and therefore

$|F_{1}(s)|\leq 2|\sigma|\sum_{i\geq 1}q^{i(1-\sigma)}\sum_{j\geq
0}\frac{P_{j}}{(1+q^{j})^{1+\sigma}}<\infty.$

Moreover, $F_{1}(0)=0$, so to the right of the fundamental strip we have the
singular expansions

$\hat{G}^{*}(s)\asymp-\frac{F^{\prime}_{1}(0)}{Ls},\quad\textup{and}\
\hat{G}^{*}(s)\asymp-p^{-2\ell\chi}\frac{\Gamma(\ell\chi)}{L}\frac{F_{1}(\ell\chi)}{s-\ell\chi},\
\textup{for}\ \ell\in\mathbb{Z}\setminus\\{0\\}.$

This leads to

$\hat{G}(n)=\frac{F^{\prime}_{1}(0)}{L}+\frac{1}{L}\sum_{\ell\in\mbox{$\mathbb{Z}$}\setminus\\{0\\}}\Gamma(\ell\chi)F_{1}(\ell\chi)(np^{2})^{-\ell\chi}+\mathcal{O}(n^{-\beta}),$

with any fixed $\beta<1$, where the constant term simplifies to

$\frac{F^{\prime}_{1}(0)}{L}=\frac{1}{L}\ln\Bigg{(}\prod_{i,k\geq
1}\frac{q^{i}+q^{k}}{q^{i}+q^{k}-pq^{i+k-1}}\Bigg{)}.$