Study of charmonium production in $b$-hadron decays and first evidence for the decay $B_s^0 \rightarrow \phi \phi \phi$

TL;DR

This paper measures charmonium production in b-hadron decays using LHCb data, and provides the first evidence for the decay B_s^0 → φφφ, including branching fractions, production ratios, and polarization observations.

Contribution

It presents the first measurement of the B_s^0 → φφφ decay and detailed ratios of charmonium states produced in b-hadron decays, expanding understanding of these processes.

Findings

First evidence for B_s^0 → φφφ decay with a branching fraction of (2.15 ± 0.54 ± 0.28 ± 0.21)×10^{-6}

Measured ratios of charmonium production in b-hadron decays, e.g., R^{χ_{c0}}_{η_c(1S)}=0.147±0.023±0.011

Observed a preferentially transverse polarization of φ mesons

Abstract

Using decays to $\phi$-meson pairs, the inclusive production of charmonium states in $b$-hadron decays is studied with $pp$ collision data corresponding to an integrated luminosity of $3.0fb^{-1}$, collected by the LHCb experiment at centre-of-mass energies of 7 and 8 TeV. Denoting by $B_C\equiv B(b\to CX)\times B(C\to\phi\phi)$ the inclusive branching fraction of a $b$ hadron to a charmonium state $C$ that decays into a pair of $\phi$ mesons, ratios $R^{C1}_{C2}\equiv B_{C1}/B_{C2}$ are determined as $R^{\chi_{c0}}_{\eta_c(1S)}=0.147\pm0.023\pm0.011$, $R^{\chi_{c1}}_{\eta_c (1S)}=0.073\pm0.016\pm0.006$, $R^{\chi_{c2}}_{\eta_c (1S)}=0.081\pm0.013\pm0.005$, $R^{\chi_{c1}}_{\chi_{c0}}=0.50\pm0.11\pm0.01$, $R^{\chi_{c2}}_{\chi_{c0}}=0.56\pm0.10\pm0.01$ and $R^{\eta_c (2S)}_{\eta_c(1S)}=0.040\pm0.011\pm0.004$. Here and below the first uncertainties are statistical and the second systematic. Upper limits at $90\%$ confidence level for the inclusive production of $X(3872)$, $X(3915)$ and $\chi_{c2}(2P)$ states are obtained as $R^{X(3872)}_{\chi_{c1}}<0.34$, $R^{X(3915)}_{\chi_{c0}}<0.12$ and $R^{\chi_{c2}(2P)}_{\chi_{c2}}<0.16$. Differential cross-sections as a function of transverse momentum are measured for the $\eta_c(1S)$ and $\chi_c$ states. The branching fraction of the decay $B_s^0\rightarrow\phi\phi\phi$ is measured for the first time, $B(B_s^0\to\phi\phi\phi)=(2.15\pm0.54\pm0.28\pm0.21_{B})\times 10^{-6}$. Here the third uncertainty is due to the branching fraction of the decay $B_s^0\to\phi\phi$, which is used for normalization. No evidence for intermediate resonances is seen. A preferentially transverse $\phi$ polarization is observed. The measurements allow the determination of the ratio of the branching fractions for the $\eta_c(1S)$ decays to $\phi\phi$ and $p\bar{p}$ as $B(\eta_c(1S)\to\phi\phi)/B(\eta_c(1S)\to p\bar{p})=1.79\pm0.14\pm0.32$.

Tables

Table 1. Table 1: Signal yields with statistical uncertainties of the fit to the spectrum of the ϕ ​ ϕ italic-ϕ italic-ϕ \phi\phi invariant mass.

Resonance	Signal yield
${\eta }_c\left(1S\right)$	6476	$\pm$	418
${\chi }_{c0}$	933	$\pm$	128
${\chi }_{c1}$	460	$\pm$	89
${\chi }_{c2}$	611	$\pm$	97
${\eta }_c\left(2S\right)$	365	$\pm$	100

Table 2. Table 2: The ratio of charmonium signal yields with respect to the η c ​ ( 1 ​ S ) subscript 𝜂 𝑐 1 𝑆 {\eta_{c}}(1S) yield and between pairs of χ c subscript 𝜒 𝑐 \chi_{c} states. The first uncertainties are statistical and the second systematic.

Resonances	Signal yield ratio
$N_{{\chi }_{c0}}/N_{{\eta }_c\left(1S\right)}$	$0.144\pm 0.022\pm 0.011$
$N_{{\chi }_{c1}}/N_{{\eta }_c\left(1S\right)}$	$0.071\pm 0.015\pm 0.006$
$N_{{\chi }_{c2}}/N_{{\eta }_c\left(1S\right)}$	$0.094\pm 0.016\pm 0.006$
$N_{{\chi }_{c1}}/N_{{\chi }_{c0}}$	$0.494\pm 0.107\pm 0.012$
$N_{{\chi }_{c2}}/N_{{\chi }_{c0}}$	$0.656\pm 0.121\pm 0.015$
$N_{{\eta }_c\left(2S\right)}/N_{{\eta }_c\left(1S\right)}$	$0.056\pm 0.016\pm 0.005$

Table 3. Table 3: Systematic uncertainties of the charmonium event yield ratios within families and with respect to the η c ​ ( 1 ​ S ) subscript 𝜂 𝑐 1 𝑆 {\eta_{c}}(1S) yield. The total uncertainty is the sum in quadrature of the individual contributions.

Systematic uncertainty	$\frac{N_{{\chi }_{c0}}}{N_{{\eta }_c\left(1S\right)}}$	$\frac{N_{{\chi }_{c1}}}{N_{{\eta }_c\left(1S\right)}}$	$\frac{N_{{\chi }_{c2}}}{N_{{\eta }_c\left(1S\right)}}$	$\frac{N_{{\chi }_{c1}}}{N_{{\chi }_{c0}}}$	$\frac{N_{{\chi }_{c2}}}{N_{{\chi }_{c0}}}$	$\frac{N_{{\eta }_c\left(2S\right)}}{N_{{\eta }_c\left(1S\right)}}$
Including other states	0.004	0.003	0.003	$0.006$	$0.008$	$0.003$
Description of detector resolution				$0.001$	$0.001$	$0.002$
Description of signal resonances	$0.002$	$0.002$	$0.001$	$0.010$	$0.002$	$0.003$
Background model	$0.010$	$0.005$	$0.005$	0.002	0.012
2D fit functions	$0.002$		$0.001$	$0.005$	$0.005$	$0.001$
Total	$0.011$	$0.006$	$0.006$	$0.012$	$0.015$	$0.005$

Table 4. Table 4: Exponential slope parameter in units of GeV − 1 superscript GeV 1 \mathrm{\,Ge\kern-1.00006ptV}^{-1} from a fit to the p T subscript 𝑝 T p_{\mathrm{T}} spectra of η c ​ ( 1 ​ S ) subscript 𝜂 𝑐 1 𝑆 {\eta_{c}}(1S) , χ c ​ 0 subscript 𝜒 𝑐 0 \chi_{{c}0} , χ c ​ 1 subscript 𝜒 𝑐 1 \chi_{{c}1} and χ c ​ 2 subscript 𝜒 𝑐 2 \chi_{{c}2} mesons.

	${\eta }_c\left(1S\right)$	${\chi }_{c0}$	${\chi }_{c1}$	${\chi }_{c2}$
$\sqrt{s}=7\mathrm{TeV}$	$0.41\pm 0.02$	$0.32\pm 0.04$	$0.31\pm 0.06$	$0.30\pm 0.05$
$\sqrt{s}=8\mathrm{TeV}$	$0.39\pm 0.02$	$0.37\pm 0.04$	$0.41\pm 0.06$	$0.33\pm 0.04$

Table 5. Table 5: Charmonium masses and natural widths in MeV MeV \,\!\mathrm{\,Me\kern-1.00006ptV} .

	Measured value	World average [14]
$M_{{\eta }_c\left(1S\right)}$	$2982.8\pm 1.0\pm 0.5$	2983.4	$\pm$	0.5
$M_{{\chi }_{c0}}$	$3413.0\pm 1.9\pm 0.6$	3414.75	$\pm$	0.31
$M_{{\chi }_{c1}}$	$3508.4\pm 1.9\pm 0.7$	3510.66	$\pm$	0.07
$M_{{\chi }_{c2}}$	$3557.3\pm 1.7\pm 0.7$	3556.20	$\pm$	0.09
$M_{{\eta }_c\left(2S\right)}$	$3636.4\pm 4.1\pm 0.7$	3639.2	$\pm$	1.2
${\Gamma }_{{\eta }_c\left(1S\right)}$	$31.4\pm 3.5\pm 2.0$	31.8	$\pm$	0.8
${\Gamma }_{{\eta }_c\left(2S\right)}$	$-$	${11.3}_{-\mathrm{\hspace{0.17em}\; 2.9}}^{+\mathrm{\hspace{0.17em}\; 3.2}}$

Table 6. Table 6: Charmonium mass differences (in MeV MeV \,\!\mathrm{\,Me\kern-1.00006ptV} ).

	Measured value	World average [14]
$M_{{\chi }_{c1}}-M_{{\chi }_{c0}}$	$95.4\pm 2.7\pm 0.1$	95.91	$\pm$	0.83
$M_{{\chi }_{c2}}-M_{{\chi }_{c0}}$	$144.3\pm 2.6\pm 0.2$	141.45	$\pm$	0.32
$M_{{\eta }_c\left(2S\right)}-M_{{\eta }_c\left(1S\right)}$	$653.5\pm 4.2\pm 0.4$	655.70	$\pm$	1.48

Equations

$$R_{C_2}^{C_1}\equiv \frac{\mathcal{B}(b\to C_{1}X)\times \mathcal{B}(C_1\to \varphi \varphi )}{\mathcal{B}(b\to C_{2}X)\times \mathcal{B}(C_2\to \varphi \varphi )}$$

$$R_{{\eta }_c(1S)}^{{\chi }_{c0}}$$

$$R_{{\eta }_c(1S)}^{{\chi }_{c1}}$$

$$R_{{\eta }_c(1S)}^{{\chi }_{c2}}$$

$$R_{{\chi }_{c0}}^{{\chi }_{c1}}$$

$$R_{{\chi }_{c0}}^{{\chi }_{c2}}$$

$$R_{{\eta }_c(1S)}^{{\eta }_c(2S)}$$

$$\frac{\mathcal{B}(b\to {\chi }_{c1}X)}{\mathcal{B}(b\to {\chi }_{c0}X)}$$

$$\frac{\mathcal{B}(b\to {\chi }_{c2}X)}{\mathcal{B}(b\to {\chi }_{c0}X)}$$

$$\frac{\mathcal{B}(b\to {\chi }_{c0}X)}{\mathcal{B}(b\to {\eta }_c(1S)X)}$$

$$\frac{\mathcal{B}(b\to {\chi }_{c1}X)}{\mathcal{B}(b\to {\eta }_c(1S)X)}$$

$$\frac{\mathcal{B}(b\to {\chi }_{c2}X)}{\mathcal{B}(b\to {\eta }_c(1S)X)}$$

$$\mathcal{B}(b\to {\chi }_{c0}X)$$

$$\mathcal{B}(b\to {\chi }_{c1}X)$$

$$\mathcal{B}(b\to {\chi }_{c2}X)$$

$$\mathcal{B}(b\to {\eta }_c(2S)X)\times \mathcal{B}({\eta }_c(2S)\to \varphi \varphi )$$

$$R_{{\chi }_{c1}}^{X(3872)}$$

$$R_{{\chi }_{c0}}^{X(3915)}$$

$$R_{{\chi }_{c2}}^{{\chi }_{c2}(2P)}$$

$$\mathcal{B}(b\to X(3872)X)\times \mathcal{B}(X(3872)\to \varphi \varphi )$$

$$\mathcal{B}(b\to X(3915)X)\times \mathcal{B}(X(3915)\to \varphi \varphi )$$

$$\mathcal{B}(b\to {\chi }_{c2}(2P)X)\times \mathcal{B}({\chi }_{c2}(2P)\to \varphi \varphi )$$

$$\frac{\mathcal{B}(B_s^0\to \varphi \varphi \varphi )}{\mathcal{B}(B_s^0\to \varphi \varphi )}=\frac{N_{B_s^0\to \varphi \varphi \varphi }}{N_{B_s^0\to \varphi \varphi }}\times \frac{{\epsilon }_{B_s^0\to \varphi \varphi }}{{\epsilon }_{B_s^0\to \varphi \varphi \varphi }}\times \frac{1}{\mathcal{B}(\varphi \to K^{+}{K}^{-})}=0.117\pm 0.030\pm \mathrm{0.015.}$$

$$\mathcal{B}(B_s^0\to \varphi \varphi \varphi )=(2.15\pm 0.54\pm 0.28\pm 0.21_{\mathcal{B}})\times 10^{-6},$$

$$R_{C_2}^{C_1}\equiv \frac{\mathcal{B}(b\to C_{1}X)\times \mathcal{B}(C_1\to \varphi \varphi )}{\mathcal{B}(b\to C_{2}X)\times \mathcal{B}(C_2\to \varphi \varphi )},$$

$$R_{{\eta }_c(1S)}^{{\chi }_{c0}}$$

$$R_{{\eta }_c(1S)}^{{\chi }_{c1}}$$

$$R_{{\eta }_c(1S)}^{{\chi }_{c2}}$$

$$R_{{\chi }_{c0}}^{{\chi }_{c1}}$$

$$R_{{\chi }_{c0}}^{{\chi }_{c2}}$$

$$R_{{\eta }_c(1S)}^{{\eta }_c(2S)}$$

$$\frac{\mathcal{B}(b\to {\chi }_{c1}X)}{\mathcal{B}(b\to {\chi }_{c0}X)}$$

$$\frac{\mathcal{B}(b\to {\chi }_{c2}X)}{\mathcal{B}(b\to {\chi }_{c0}X)}$$

$$\mathcal{B}(b\to {\chi }_{c0}X)$$

$$\mathcal{B}(b\to {\chi }_{c1}X)$$

$$\mathcal{B}(b\to {\chi }_{c2}X)$$

$$\mathcal{B}(b\to {\eta }_c(2S)X)\times \mathcal{B}({\eta }_c(2S)\to \varphi \varphi )$$

$$R_{{\chi }_{c1}}^{X(3872)}<0.39(0.34),R_{{\chi }_{c0}}^{X(3915)}<0.14(0.12)\mathrm{and}{R}_{{\chi }_{c2}}^{{\chi }_{c2}(2P)}<0.20(0.16),$$

$$\mathcal{B}(b\to X(3872)X)\times \mathcal{B}(X(3872)\to \varphi \varphi )$$

$$\mathcal{B}(b\to X(3915)X)\times \mathcal{B}(X(3915)\to \varphi \varphi )$$

$$\mathcal{B}(b\to {\chi }_{c2}(2P)X)\times \mathcal{B}({\chi }_{c2}(2P)\to \varphi \varphi )$$

$$\mathcal{B}(B_s^0\to \varphi \varphi \varphi )=(2.15\pm 0.54\pm 0.28\pm 0.21_{\mathcal{B}})\times 10^{-6}.$$

Full text

EUROPEAN ORGANIZATION FOR NUCLEAR RESEARCH (CERN)

​​​

CERN-EP-2017-099

LHCb-PAPER-2017-007

Study of charmonium production in $b$-hadron decays and first evidence for the decay ${{B}^{0}_{s}}\!\rightarrow\phi\phi\phi$

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

Using decays to $\phi$-meson pairs, the inclusive production of charmonium states in $b$-hadron decays is studied with $pp$ collision data corresponding to an integrated luminosity of $3.0\mbox{\,fb}^{-1}$, collected by the LHCb experiment at centre-of-mass energies of 7 and 8 TeV. Denoting by ${\mathcal{B}}_{C}\equiv{\mathcal{B}}({b}\!\rightarrow CX)\times{\mathcal{B}}(C\!\rightarrow\phi\phi)$ the inclusive branching fraction of a ${b}$ hadron to a charmonium state $C$ that decays into a pair of $\phi$ mesons, ratios $R^{C_{1}}_{C_{2}}\equiv{\mathcal{B}}_{C_{1}}/{\mathcal{B}}_{C_{2}}$ are determined as $R^{{\chi_{{c}0}}}_{{\eta_{c}}(1S)}=0.147\pm 0.023\pm 0.011$, $R^{{\chi_{{c}1}}}_{{\eta_{c}}(1S)}=0.073\pm 0.016\pm 0.006$, $R^{{\chi_{{c}2}}}_{{\eta_{c}}(1S)}=0.081\pm 0.013\pm 0.005$, $R^{{\chi_{{c}1}}}_{{\chi_{{c}0}}}=0.50\pm 0.11\pm 0.01$, $R^{{\chi_{{c}2}}}_{{\chi_{{c}0}}}=0.56\pm 0.10\pm 0.01$ and $R^{{\eta_{c}}(2S)}_{{\eta_{c}}(1S)}=0.040\pm 0.011\pm 0.004$. Here and below the first uncertainties are statistical and the second systematic. Upper limits at 90% confidence level for the inclusive production of $X(3872)$, $X(3915)$ and ${\chi_{{c}2}}(2P)$ states are obtained as $R^{X(3872)}_{{\chi_{{c}1}}}<0.34$, $R^{X(3915)}_{{\chi_{{c}0}}}<0.12$ and $R^{{\chi_{{c}2}}(2P)}_{{\chi_{{c}2}}}<0.16$. Differential cross-sections as a function of transverse momentum are measured for the ${\eta_{c}}(1S)$ and $\chi_{c}$ states. The branching fraction of the decay ${{B}^{0}_{s}}\!\rightarrow\phi\phi\phi$ is measured for the first time, ${\mathcal{B}}({{B}^{0}_{s}}\!\rightarrow\phi\phi\phi)=(2.15\pm 0.54\pm 0.28\pm 0.21_{{\mathcal{B}}})\times 10^{-6}$. Here the third uncertainty is due to the branching fraction of the decay ${{B}^{0}_{s}}\!\rightarrow\phi\phi$, which is used for normalization. No evidence for intermediate resonances is seen. A preferentially transverse $\phi$ polarization is observed. The measurements allow the determination of the ratio of the branching fractions for the ${\eta_{c}}(1S)$ decays to $\phi\phi$ and ${p}{\overline{p}}$ as ${\mathcal{B}}({\eta_{c}}(1S)\!\rightarrow\phi\phi)/{\mathcal{B}}({\eta_{c}}(1S)\!\rightarrow{p}{\overline{p}})=1.79\pm 0.14\pm 0.32$.

Published in Eur. Phys. J. C77 (2017) 609

© CERN on behalf of the LHCb collaboration, licence CC-BY-4.0.

1 Introduction

The production of the $J^{PC}=1^{--}$ charmonium states has been extensively studied using decays to clean dilepton final states. Other states such as those from the $\chi_{c}$ family can be accessed via the radiative transition to ${J\mskip-3.0mu/\mskip-2.0mu\psi\mskip 2.0mu}$. Studies of the production of the non-$1^{--}$ charmonium states can be performed by reconstructing their decays to fully hadronic final states [1]. This paper reports a measurement of the inclusive production rates of the $\eta_{c}$ and $\chi_{c}$ states in $b$-hadron decays, ${b}\!\rightarrow{\eta_{c}}X$ and ${b}\!\rightarrow\chi_{c}X$, using charmonia decays to a pair of $\phi$ mesons. In addition, the first evidence for the decay ${{B}^{0}_{s}}\!\rightarrow\phi\phi\phi$ is reported.

Results on inclusive charmonium production in $b$-hadron decays are available from $e^{+}e^{-}$ experiments operating at centre-of-mass energies around the ${\mathchar 28935\relax{(4S)}}$ and ${\mathchar 28935\relax{(5S)}}$ resonances, studying mixtures of ${{B}^{+}}$ and ${B}^{0}$ mesons111The inclusion of charge-conjugate states is implied throughout. (light mixture) or ${{B}^{+}}$, ${B}^{0}$ and ${B}^{0}_{s}$ mesons, respectively. Mixtures of all $b$-hadrons (${{B}^{+}}$, ${B}^{0}$, ${B}^{0}_{s}$, ${B}_{c}^{+}$ and $b$-baryons) have been studied at LEP, the Tevatron and the LHC. The world average values for charmonium branching fractions from the light mixture are dominated by results from the CLEO [2, 3], Belle [4] and BaBar [5] collaborations. For the ${J\mskip-3.0mu/\mskip-2.0mu\psi\mskip 2.0mu}$, $\psi{(2S)}$ and $\chi_{{c}1}$ states the measured branching fractions are consistent within uncertainties. The new Belle result for the ${b}\!\rightarrow{\chi_{{c}2}}X$ branching fraction [4], which supersedes the previous measurement [6], is below the BaBar result [5] by more than 2.5 standard deviations, while the CLEO collaboration does not observe a statistically significant ${b}\!\rightarrow{\chi_{{c}2}}X$ signal [3]. An upper limit on the inclusive production rate of ${\eta_{c}}(1S)$ mesons in the light mixture, ${\mathcal{B}}({B}\rightarrow{\eta_{c}}(1S)X)<9\times 10^{-3}$ at $90\%$ confidence level (CL), was reported by CLEO [7].

The branching fractions of $b$-hadron decays to final states including a ${J\mskip-3.0mu/\mskip-2.0mu\psi\mskip 2.0mu}$ or $\psi{(2S)}$ charmonium state, where all $b$-hadron species are involved, are known with uncertainties of around 10%, with the world averages dominated by the measurements performed at LEP [8, 9, 10]. The ratio of ${b}\!\rightarrow{\psi{(2S)}}X$ and ${b}\!\rightarrow{{J\mskip-3.0mu/\mskip-2.0mu\psi\mskip 2.0mu}}X$ yields has been measured at the LHC by the LHCb, CMS and ATLAS collaborations with a precision of around 5% [11, 12, 13]. The only available results for the $\chi_{c}$ family are the $\chi_{{c}1}$ inclusive production rates in $b$-hadron decays measured by the DELPHI and L3 collaborations [8, 9], with an average value of ${\mathcal{B}}({b}\rightarrow{\chi_{{c}1}}X)=(14\pm 4)\times 10^{-3}$ [14]. Recently, LHCb measured the ${\eta_{c}}(1S)$ production rate, ${\mathcal{B}}({b}\rightarrow{\eta_{c}}(1S)X)=(4.88\pm 0.64\pm 0.29\pm 0.67_{{\mathcal{B}}})\times 10^{-3}$, where the third uncertainty is due to uncertainties on the ${J\mskip-3.0mu/\mskip-2.0mu\psi\mskip 2.0mu}$ inclusive branching fraction from $b$-hadron decays and the branching fractions of the decays of ${J\mskip-3.0mu/\mskip-2.0mu\psi\mskip 2.0mu}$ and ${\eta_{c}}(1S)$ to the $p$ $\overline{p}$ final state [15].

While experimentally the reconstruction of charmonia from $b$-hadron decays allows an efficient control of combinatorial background with respect to charmonium candidates produced in the $p$ $p$ collision vertex via hadroproduction or in the decays of heavier resonances (prompt charmonium), inclusive $b$-hadron decays to charmonia are theoretically less clean. Since a description of the strong interaction dynamics in $b$-hadron inclusive decays improves with respect to exclusive decays due to consideration of more final states, and the formation of charmonium proceeds through a short-distance process, a factorization of a $c\bar{c}$ pair production and its hadronization in a given charmonium state becomes a reasonable assumption [16]. The relative inclusive production of $\chi_{c}$ states in $b$-hadron decays provides a clean test of charmonia production models. For example, the colour evaporation model predicts a ${\chi_{{c}2}}/{\chi_{{c}1}}$ production ratio of $5/3$ [17], while the perturbative QCD-based computation predicts that the V-A current, which is responsible for the $b$ decays, forbids the $\chi_{{c}2}$ and $\chi_{{c}0}$ production at leading order. In the non-relativistic QCD (NRQCD) framework [18, 19, 20], the colour-octet contributions have to be included, predicting the rates to be proportional to $(2J+1)$ for the ${\chi_{c}}_{J}$ states. The NRQCD framework can be applied to both prompt charmonium production and secondary production from $b$-hadron decays and the comparison between these two production mechanisms can provide a valuable test of this theoretical framework.

In this paper we report the first measurements of inclusive $\chi_{c}$ and ${\eta_{c}}(2S)$ production rates in $b$-hadron decays using charmonium decays to hadronic final states in the high-multiplicity environment of a hadron collider. Experimentally, charmonium candidates from $b$-hadron decays are distinguished from prompt charmonia by exploiting the $b$-hadron decay time and reconstructing a $b$-hadron (and charmonium) decay vertex well separated from the primary vertex where the $b$-hadron candidate was produced. The charmonium states are reconstructed via their decays to a $\phi\phi$ final state. The ${\eta_{c}}(1S)$ production followed by the decay ${\eta_{c}}(1S)\rightarrow\phi\phi$ is used for normalization, so that systematic uncertainties partially cancel in the ratios. As a by-product of the production rate measurements, the masses of the ${\eta_{c}}(1S)$, $\chi_{{c}0}$, $\chi_{{c}1}$, $\chi_{{c}2}$ and ${\eta_{c}}(2S)$ charmonium states and the natural width of the ${\eta_{c}}(1S)$ meson are determined.

The ${B}^{0}_{s}$ decay to the $\phi\phi$ final state has been observed by the CDF collaboration [21] and recently precisely measured by the LHCb collaboration [22], where it was also used to search for $C\!P$-violating asymmetries [23]. In the Standard Model (SM) the amplitude for the decay ${{B}^{0}_{s}}\rightarrow\phi\phi$ is dominated by a loop diagram. Experimental verification of the partial width, polarization amplitudes and triple-product asymmetries of the ${{B}^{0}_{s}}\rightarrow\phi\phi$ decay probes the QCD contribution to the weak processes described by nonfactorizable penguin diagrams [24, 25], and contributions from particles beyond the SM to the penguin loops [26, 27, 28, 29, 30]. A three-body ${{B}^{0}_{s}}\rightarrow\phi\phi\phi$ decay leads to a final state with six strange quarks. In the SM it is described by the penguin diagram of the ${{B}^{0}_{s}}\rightarrow\phi\phi$ decay with the creation of an additional ${s}{\overline{s}}$ quark pair. The ${{B}^{0}_{s}}\rightarrow\phi\phi\phi$ decay can also receive contributions from an intermediate charmonium state decaying to a $\phi\phi$ state. Here we report first evidence for the ${{B}^{0}_{s}}\rightarrow\phi\phi\phi$ decay and study its resonance structure. The branching fraction of this decay is determined relative to the branching fraction ${\mathcal{B}}({{B}^{0}_{s}}\rightarrow\phi\phi)$ [22]. To cross-check the technique exploited in this paper, the value of ${\mathcal{B}}({{B}^{0}_{s}}\rightarrow\phi\phi)$ is also determined relative to the ${\eta_{c}}(1S)$ production rate. Finally, the ratio of the branching fractions for the decays ${\eta_{c}}(1S)\rightarrow\phi\phi$ and ${\eta_{c}}(1S)\rightarrow{p}{\overline{p}}$ is determined, using additional external information.

The LHCb detector and data sample used for the analysis are presented in Sect. 2. Section 3 explains the selection details and the signal extraction technique. Inclusive production of charmonium states in $b$-hadron decays is discussed in Sect. 4. In Sect. 5 measurements of the ${\eta_{c}}(1S)$ mass and natural width are described. First evidence for the ${{B}^{0}_{s}}\rightarrow\phi\phi\phi$ decay is reported in Sect. 6. The main results of the paper are summarized in Sect. 7.

2 LHCb detector and data sample

The LHCb detector [31, 32] is a single-arm forward spectrometer covering the pseudorapidity range $2<\eta<5$, designed for the study of particles containing $b$ or $c$ quarks. The detector includes a high-precision tracking system consisting of a silicon-strip vertex detector surrounding the $pp$ interaction region, a large-area silicon-strip detector located upstream of a dipole magnet with a bending power of about $4{\mathrm{\,Tm}}$, and three stations of silicon-strip detectors and straw drift tubes placed downstream of the magnet. The tracking system provides a measurement of momentum, $p$, of charged particles with a relative uncertainty that varies from 0.5% at low momentum222Natural units are used throughout the paper. to 1.0% at 200$\mathrm{\,Ge\kern-1.00006ptV}$. The minimum distance of a track to a primary vertex (PV), the impact parameter (IP), is measured with a resolution of (15+29/\mbox{p_{\mathrm{T}}}){\,\upmu\mathrm{m}}, where $p_{\mathrm{T}}$ is the component of the momentum transverse to the beam, in $\mathrm{\,Ge\kern-1.00006ptV}$. Different types of charged hadrons are distinguished using information from two ring-imaging Cherenkov detectors. Photons, electrons and hadrons are identified by a calorimeter system consisting of scintillating-pad and preshower detectors, an electromagnetic calorimeter and a hadronic calorimeter. Muons are identified by a system composed of alternating layers of iron and multiwire proportional chambers. 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.

The analysis is based on $p$ $p$ collision data recorded by the LHCb experiment at a centre-of-mass energy $\sqrt{s}=7\mathrm{\,Te\kern-1.00006ptV}$, corresponding to an integrated luminosity of $1.0\mbox{\,fb}^{-1}$, and at $\sqrt{s}=8\mathrm{\,Te\kern-1.00006ptV}$, corresponding to an integrated luminosity of $2.0\mbox{\,fb}^{-1}$. Events enriched in signal decays are selected by the hardware trigger, based on the presence of a single deposit of high transverse energy in the calorimeter. The subsequent software trigger selects events with displaced vertices formed by charged particles having a good track-fit quality, transverse momentum larger than 0.5$\mathrm{\,Ge\kern-1.00006ptV}$, and that are incompatible with originating from any PV [23]. Charged kaon candidates are identified using the information from the Cherenkov and tracking detectors. Two oppositely charged kaon candidates having an invariant mass within $\pm 11\mathrm{\,Me\kern-1.00006ptV}$ of the known mass of the $\phi$ meson are required to form a good quality vertex.

Precise mass measurements require a momentum-scale calibration. The procedure [33] uses ${{J\mskip-3.0mu/\mskip-2.0mu\psi\mskip 2.0mu}}\rightarrow{\mu^{+}}{\mu^{-}}$ decays to cross-calibrate a relative momentum scale between different data-taking periods. The absolute scale is determined using ${{{B}^{+}}}\rightarrow{{J\mskip-3.0mu/\mskip-2.0mu\psi\mskip 2.0mu}}{{K}^{+}}$ decays with known particle masses as input [14]. The final calibration is checked with a variety of fully reconstructed quarkonium, ${{B}^{+}}$ and ${K}^{0}_{\mathrm{\scriptscriptstyle S}}$ decays. No residual bias is observed within the experimental resolution.

In the simulation, $pp$ collisions are generated using Pythia [34, 35] with a specific LHCb configuration [36]. Decays of hadronic particles are described by EvtGen [37], in which final-state radiation is generated using Photos [38]. The interaction of the generated particles with the detector material and its response are implemented using the Geant4 toolkit [39, *Agostinelli:2002hh] as described in Ref. [41]. Simulated samples of ${\eta_{c}}$ and $\chi_{c}$ mesons decaying to the $\phi\phi$ final state, and ${B}^{0}_{s}$ decaying to two and three $\phi$ mesons, are used to estimate efficiency ratios and to evaluate systematic uncertainties. Charmonium and ${{B}^{0}_{s}}\rightarrow\phi\phi\phi$ decays are generated with uniform phase-space density, while ${{B}^{0}_{s}}\rightarrow\phi\phi$ decays are generated according to the amplitudes from Ref. [21].

3 Selection criteria and signal extraction

The signal selection is largely performed at the trigger level. The offline analysis selects combinations of two or three $\phi$ candidates that are required to form a good-quality common vertex displaced from the primary vertex. A good separation between the two vertices ($\chi^{2}>100$) is required, reducing the contribution from charmonia produced directly at the primary vertex to below $1\%$. Pairs of $\phi$ mesons originating from different $b$-hadrons produced in the same beam crossing event are suppressed by the requirement of a good-quality common vertex. Detector acceptance and selection requirements, and in particular the requirement of the kaon $p_{\mathrm{T}}$ to exceed 0.5$\mathrm{\,Ge\kern-1.00006ptV}$, significantly suppress $\phi\phi$ combinations with $p_{\mathrm{T}}$ below $4\mathrm{\,Ge\kern-1.00006ptV}$.

Two-dimensional (2D) or three-dimensional (3D) extended unbinned maximum likelihood fits, corresponding to the two or three randomly ordered $K^{+}K^{-}$ combinations, are performed in bins of the invariant mass of the four-kaon and six-kaon combinations, denoted as $2(K^{+}K^{-})$ or $3(K^{+}K^{-})$, respectively, to determine the numbers of $\phi\phi$ and $\phi\phi\phi$ combinations. The 2D fit accounts for the signal, $\phi\phi$, and background, $\phi\,(K^{+}K^{-})$ and $2(K^{+}K^{-})$, components, while the 3D fit accounts for the signal, $\phi\phi\phi$, and background, $\phi\phi\,(K^{+}K^{-})$, $\phi\,2(K^{+}K^{-})$ and $3(K^{+}K^{-})$, components. A $\phi$ signal is described by the convolution of a Breit-Wigner function and a sum of two Gaussian functions with a common mean. The ratio of the two Gaussian widths and the fraction of the narrower Gaussian are taken from simulation. The contribution from combinatorial background, due to ${{K}^{+}}{{K}^{-}}$ pairs not originating from the decay of a $\phi$ meson, is assumed to be flat. In addition, a threshold function to account for the available phase-space in the ${{K}^{+}}{{K}^{-}}$ system is introduced for both signal and combinatorial background. While no visible contribution from the $f_{0}(980)$ resonance decaying into a $K^{+}K^{-}$ pair is observed in the $2(K^{+}K^{-})$ or $3(K^{+}K^{-})$ combinations, a potential effect due to contributions from such decays is considered as a source of systematic uncertainty.

Figures 1 and 2 show the results of the 2D fits to the $2(K^{+}K^{-})$ invariant mass distributions along with the projections to the ${{K}^{+}}{{K}^{-}}$ invariant mass axes in the ${\eta_{c}}(1S)$ and ${B}^{0}_{s}$ signal regions, $2.91-3.06\mathrm{\,Ge\kern-1.00006ptV}$ and $5.30-5.43\mathrm{\,Ge\kern-1.00006ptV}$. Figure 3 shows the projections to the ${{K}^{+}}{{K}^{-}}$ invariant mass axes of the 3D fit to the $3(K^{+}K^{-})$ invariant mass distribution in the ${B}^{0}_{s}$ signal region. While using the known value for the natural width of the $\phi$ resonance [14], the $\phi$ mass and the remaining resolution parameter are determined from the fits in the enlarged signal $\phi\phi$ and $\phi\phi\phi$ invariant mass regions. In the 2D and 3D fits in the bins of $\phi\phi$ and $\phi\phi\phi$ invariant mass the $\phi$ mass and the resolution parameter are then fixed to the values determined in the enlarged signal regions.

Unless they are extracted from the 2D or 3D fits, throughout the paper the error bars shown in the histograms are determined as follows: the upper (lower) error bar assigned to a given bin content $N$ is defined by the expectation value $\lambda$ of the Poissonian distribution giving 16% probability to observe the number of events $N$ or less (more). When obtained from the 2D or 3D fits the histogram bin contents are constrained to be positive, with error bars defined by the range in the allowed region where the best fit negative-log likelihood value is within half a unit from the global minimum.

In the following, production ratios are determined from the signal yields obtained from the fits of the $\phi\phi$ or $\phi\phi\phi$ invariant mass spectra. The relative efficiencies are taken into account to determine the ratio of the branching fractions of the corresponding decays. Signal yields corresponding to the data samples accumulated at $\sqrt{s}=7$ and $8\mathrm{\,Te\kern-1.00006ptV}$ are found to be compatible. Unless otherwise specified, the results below are based on the analysis of the combined data sample.

4 Charmonium production in decays to $\boldsymbol{\phi\phi}$

4.1 Charmonium yields

Figure 4 shows the invariant mass spectrum of the $\phi\phi$ combinations, where the content of each bin is a result of a 2D fit to the two ${{K}^{+}}{{K}^{-}}$ invariant-mass combinations.

A binned $\chi^{2}$ fit to the spectrum is used to determine the contributions from the ${\eta_{c}}(1S)$ and ${\eta_{c}}(2S)$ mesons, and the $\chi_{{c}0}$, $\chi_{{c}1}$ and $\chi_{{c}2}$ mesons. The charmonium-like states $X(3872)$, $X(3915)$ and ${\chi_{{c}2}}(2P)$ with masses and natural widths from Ref. [14] are taken into account in alternative fits in order to evaluate systematic uncertainties, as well as to obtain upper limits on the inclusive production of these states in $b$-hadron decays. Each resonance is described by the convolution of a relativistic Breit-Wigner function and a sum of two Gaussian functions with a common mean. The combinatorial background, i.e. contributions due to random combinations of two true $\phi$ mesons, is described by the product of a first-order polynomial with an exponential function and a threshold factor. The natural width of the ${\eta_{c}}(1S)$ state is a free parameter in the fit, while the natural widths of the ${\eta_{c}}(2S)$ and the $\chi_{c}$ states, which have lower signal yields, are fixed to their known values [14]. Possible variations of the ${\eta_{c}}(2S)$ production rate depending on its natural width, $\Gamma_{{\eta_{c}}(2S)}$, are explored by providing the result as a function of the assumed natural width. The ratio of the two Gaussian widths and the fraction of the narrow Gaussian are fixed from the simulation. The mass resolution for different charmonium resonances is scaled according to the energy release, so that a single free parameter in the $\phi\phi$ invariant mass fit accounts for the detector resolution. This scaling of the mass resolution for different charmonium states has been validated using simulation.

The signal yields are given in Table 1. The ratios of the resonance yields from the fit are given in Table 2, both for the ratios with respect to the ${\eta_{c}}(1S)$ yield and between pairs of $\chi_{c}$ states; the systematic uncertainties are discussed below. The statistical significance for the $N_{{\eta_{c}}(2S)}$ signal is estimated from the $\chi^{2}$-profile to be 3.7 standard deviations.

Systematic uncertainties in the ratios of the charmonium yields are estimated by considering potential contributions from other states, from imperfect modelling of detector resolution, signal resonances and background, and from the 2D fit technique. In order to evaluate the systematic uncertainty related to potential contributions from other states, signal shapes for the $X(3872)$, $X(3915)$, and ${\chi_{{c}2}}(2P)$ states are included in the fit. Systematic uncertainties related to detector resolution are estimated by fixing the ${\eta_{c}}(1S)$ mass resolution to the value determined from the simulation. In addition, systematic uncertainties associated to the impact of the detector resolution on the signal shapes are estimated by comparing the nominal fit results to those obtained using a single instead of a double Gaussian function. An uncertainty associated with the description of the mass resolution of the $\phi$ meson is estimated by fixing the resolution in the 2D fits to the value determined from simulation. The uncertainty associated with the description using the relativistic Breit-Wigner line shape [42] is estimated by varying the radial parameter of the Blatt-Weisskopf barrier factor [43] between $0.5$ and $3\,\mathrm{\,Ge\kern-1.00006ptV}^{-1}$. In order to estimate the uncertainty related to the natural width of the ${\eta_{c}}(2S)$ meson, the value of $\Gamma_{{\eta_{c}}(2S)}$ is varied within the uncertainties of the world average [14]. The uncertainty in the description of the $\chi_{c}$ signal peaks is estimated by fixing the $\chi_{c}$ masses to their known values. A reduced fit range, covering only the $\chi_{c}$ and ${\eta_{c}}(2S)$ region ($3.15-3.95\mathrm{\,Ge\kern-1.00006ptV}$), is used to estimate a systematic uncertainty associated to the choice of the fit range. An alternative background parametrization using a parabolic instead of a linear function is used to estimate the systematic uncertainty due to the choice of the background parametrisation. A systematic uncertainty associated to the background parametrization in the 2D fits is estimated by adding slope parameters to the description of the non-$\phi$ ${{K}^{+}}{{K}^{-}}$ combinations in the $\phi{{K}^{+}}{{K}^{-}}$ and the $2\times({{K}^{+}}{{K}^{-}})$ components. The effect of a potential contribution from the $f_{0}(980)$ state in the 2D fits is estimated by including the $f_{0}(980)$ contribution following the analysis described in Ref. [44], and varying the $f_{0}(980)$ mass and natural width within the uncertainties quoted in Ref. [14]. Contributions from multiple candidates are below 2% per event and are not considered in the evaluation of systematic uncertainties. The uncertainty related to the momentum-scale calibration is negligible.

The total systematic uncertainty is obtained as the quadratic sum of the individual systematic contributions. The systematic uncertainties are shown in Table 3.

The description of the background and the potential contributions from other resonances dominate the total systematic uncertainties. In the yield ratios the systematic uncertainty is smaller than or comparable to the statistical uncertainty.

Complementary cross-checks are performed by comparing the distributions of kinematic variables in simulation and data. The stability of the results is checked by using an alternative binning of the $\phi\phi$ invariant mass distribution (shifted by half a bin) and by using the weighted signal events from the sPlot [45] instead of the nominal 2D fit technique. The same check is performed for the determination of the masses and widths of the charmonium states. In all cases no significant changes are observed and no additional contributions to the systematic uncertainties are assigned.

4.2 Production of $\boldsymbol{{\eta_{c}}}$ and $\boldsymbol{\chi_{c}}$ in $b$-hadron decays

The production ratios of charmonium $C$ with respect to the ${\eta_{c}}(1S)$ yield and between pairs of $\chi_{c}$ states

[TABLE]

are determined to be

[TABLE]

where, here and throughout, the first uncertainties are statistical and the second are systematic. The dominant contributions to the systematic uncertainty on the relative $\chi_{c}$ production rates arise from accounting for possible other resonances and from uncertainties on the $\chi_{c}$ masses [14]. The systematic uncertainties are smaller than the statistical uncertainties, so that the overall uncertainty on the measurements will be reduced with a larger dataset. The systematic uncertainty on the relative production rate of ${\eta_{c}}(2S)$ mesons is dominated by possible contributions from other resonances and variation of the ${\eta_{c}}(2S)$ natural width.

In the following, the ${\eta_{c}}(1S)$ production rate in $b$-hadron decays and the branching fractions of the charmonium decays to $\phi\phi$ are used to obtain single ratios of charmonium production rates in $b$-hadron decays and the branching fractions of inclusive $b$-hadron transitions to charmonium states. The ${\eta_{c}}(1S)$ inclusive production rate in $b$-hadron decays was measured by LHCb as ${\mathcal{B}}(b\rightarrow{\eta_{c}}(1S)X)=(4.88\pm 0.97)\times 10^{-3}$ [15] using decays to $p$ $\overline{p}$. Branching fractions of the charmonia decays to $\phi\phi$ from Ref. [14] are used. However, Ref. [14] indicates a possible discrepancy for the ${\mathcal{B}}({\eta_{c}}(1S)\rightarrow\phi\phi)$ value when comparing a direct determination and a fit including all available measurements. Therefore, an average of the results from Belle [46] and BaBar [47] using ${{B}^{+}}$ decays to $\phi\phi{{K}^{+}}$, ${\mathcal{B}}({\eta_{c}}(1S)\rightarrow\phi\phi)=(3.21\pm 0.72)\times 10^{-3}$, is used below. The uncertainty of this average dominates the majority of the following results in this section, and an improved knowledge of the ${\mathcal{B}}({\eta_{c}}(1S)\rightarrow\phi\phi)$ is critical to reduce the uncertainties of the derived results. The values ${\mathcal{B}}({\chi_{{c}0}}\rightarrow\phi\phi)=(7.7\pm 0.7)\times 10^{-4}$, ${\mathcal{B}}({\chi_{{c}1}}\rightarrow\phi\phi)=(4.2\pm 0.5)\times 10^{-4}$, and ${\mathcal{B}}({\chi_{{c}2}}\rightarrow\phi\phi)=(1.12\pm 0.10)\times 10^{-3}$, are used for the $\chi_{c}$ decays [14]. The relative branching fractions of $b$-hadron inclusive decays into $\chi_{c}$ states are then derived to be

[TABLE]

where the third uncertainty is due to those on the branching fractions ${\mathcal{B}}(\chi_{c}\rightarrow\phi\phi)$. The result for the relative $\chi_{{c}1}$ and $\chi_{{c}2}$ production in inclusive $b$-hadron decays is close to the values measured in ${B}^{0}$ and ${{B}^{+}}$ production [14].

The branching fractions of $b$-hadron decays into $\chi_{c}$ states relative to the ${\eta_{c}}(1S)$ meson are

[TABLE]

where the dominating uncertainty is due to the uncertainty of the branching fractions ${\mathcal{B}}({\eta_{c}}(1S)\rightarrow\phi\phi)$ and ${\mathcal{B}}(\chi_{c}\rightarrow\phi\phi)$. The absolute branching fractions are determined to be

[TABLE]

where the third uncertainty is due to the uncertainties on the branching fractions of the $b$-hadron decays to the ${\eta_{c}}(1S)$ meson, ${\mathcal{B}}({b}\rightarrow{\eta_{c}}(1S)X)$, and of ${\eta_{c}}(1S)$ and $\chi_{c}$ decays to $\phi\phi$. The branching fraction of $b$-hadron decays into $\chi_{{c}0}$ is measured for the first time, and is found to be larger than the values predicted in Ref. [18].

Throughout the paper comparisons of the results on the production of charmonium states to theory predictions neglect the fact that the measured branching fractions also contain decays via intermediate higher-mass charmonium resonances, whereas theory calculations consider only direct $b$-hadron transitions to the charmonium state considered. Among the contributions that can be quantified the most sizeable comes from the $\psi{(2S)}$ state decaying to the $\chi_{c}$ states. With the branching fraction ${\mathcal{B}}({b}\!\rightarrow{\psi{(2S)}}X)$ recently measured [11] by LHCb the branching fractions ${\mathcal{B}}({b}\!\rightarrow\chi_{c}X)$, measured in this paper, are influenced by about $10\%$. The branching fractions ${\mathcal{B}}({b}\!\rightarrow{\chi_{{c}0}}X)$ and ${\mathcal{B}}({b}\!\rightarrow{\chi_{{c}2}}X)$ remain different from the predictions in Ref. [18].

The branching fraction measurement for $b$-hadron decays into $\chi_{{c}1}$ is most precise in mixtures of ${B}^{0}$, ${{B}^{+}}$, ${B}^{0}_{s}$, ${B}_{c}^{+}$ and $b$ baryons. The central value is lower than the central values measured by the DELPHI [8] and L3 [9] experiments at LEP, $(11.3^{+5.8}_{-5.0}\pm 0.4)\times 10^{-3}$ and $(19\pm 7\pm 1)\times 10^{-3}$, respectively. For the measurements with different $b$-hadron content, the LHCb result is consistent with measurements by CLEO [2], Belle [4], and BaBar [5]. Finally, the LHCb result for the inclusive $b$-hadron decays into $\chi_{{c}1}$ is consistent with the prediction in Ref. [18].

The branching fraction of $b$-hadron decays into $\chi_{{c}2}$ is measured for the first time with a mixture of ${B}^{0}$, ${{B}^{+}}$, ${B}^{0}_{s}$, ${B}_{c}^{+}$ and $b$-baryons. The result is consistent with the world average [14] measured with the ${B}^{0}$ and ${{B}^{+}}$ mixture, and with individual results from CLEO [3], Belle [4] and BaBar [5]. The value obtained is below the range predicted in Ref. [18].

A deviation of the ${\eta_{c}}(2S)$ natural width from the world average value [14] would affect the measured ratio of ${\eta_{c}}(2S)$ and ${\eta_{c}}(1S)$ production rates in $b$-hadron inclusive decays, as shown in Fig. 5.

The decay ${\eta_{c}}(2S)\rightarrow\phi\phi$ has not been observed so far. Hence the product of the branching fraction of $b$-hadron decays to ${\eta_{c}}(2S)$ and the branching fraction of the ${\eta_{c}}(2S)\rightarrow\phi\phi$ decay mode is determined as

[TABLE]

where the systematic uncertainty is dominated by the uncertainty on the ${\eta_{c}}(1S)$ production rate in $b$-hadron decays. This is the first evidence for ${\eta_{c}}(2S)$ production in $b$-hadron decays, and for the decay of the ${\eta_{c}}(2S)$ meson into a pair of $\phi$ mesons.

4.3 Transverse momentum dependence of the differential cross-sections for $\boldsymbol{{\eta_{c}}(1S)}$ and $\boldsymbol{\chi_{c}}$ production

The shapes of the differential production cross-sections as a function of transverse momentum are studied in the LHCb acceptance ($2<\eta<5$) and for 3<\mbox{p_{\mathrm{T}}}<17\mathrm{\,Ge\kern-1.00006ptV} and 2<\mbox{p_{\mathrm{T}}}<19\mathrm{\,Ge\kern-1.00006ptV} for the ${\eta_{c}}(1S)$ and $\chi_{c}$ states, respectively. Each differential production cross-section is normalized to the production cross-section integrated over the studied $p_{\mathrm{T}}$ region. Figure 6 shows the normalized differential cross-sections of ${\eta_{c}}(1S)$, $\chi_{{c}0}$, $\chi_{{c}1}$ and $\chi_{{c}2}$ production at $\sqrt{s}=7$ and $8\mathrm{\,Te\kern-1.00006ptV}$. An exponential function proportional to \exp(-\alpha\,\mbox{p_{\mathrm{T}}}) is fitted to the integral of the each bin of the distributions. No significant difference is observed between the $\sqrt{s}=7\mathrm{\,Te\kern-1.00006ptV}$ and $8\mathrm{\,Te\kern-1.00006ptV}$ data. The results for the slope parameters $\alpha$ are given in Table 4. For $\chi_{{c}1}$ and $\chi_{{c}2}$ production in $b$-hadron decays these results extend the ATLAS studies [48] in $p_{\mathrm{T}}$ and rapidity.

4.4 Search for production of $\boldsymbol{X(3872)}$, $\boldsymbol{X(3915)}$ and $\boldsymbol{{\chi_{{c}2}}(2P)}$

The observation of the $X(3915)$ and ${\chi_{{c}2}}(2P)$ states in $b$-hadron decays or the $X(3872)$ decaying to a pair of $\phi$ mesons would provide interesting information on the properties of these states. The invariant mass spectrum of $\phi\phi$ combinations in Fig. 4 shows no evidence for a signal from the $X(3872)$, $X(3915)$, or ${\chi_{{c}2}}(2P)$ states. Bayesian upper limits assuming a uniform prior in the event yields are obtained on the branching fractions relative to those involving decays to the states with similar quantum numbers. For the states with similar quantum numbers, in the efficiency ratios systematic uncertainties largely cancel. Using the efficiency ratios from the simulation, the upper limits at $95\%$ ($90\%$) CL on the ratios of inclusive branching fractions are

[TABLE]

Using the measured production rates of the $\chi_{c}$ states in $b$-hadron decays and branching fractions for the $\chi_{c}$ decays to the $\phi\phi$ final state [14], the upper limits at $95\%$ ($90\%$) CL on the production rates of the $X(3872)$, $X(3915)$, and ${\chi_{{c}2}}(2P)$ states in $b$-hadron decays are

[TABLE]

5 Masses and natural widths of charmonium states

The majority of the ${\eta_{c}}(1S)$ mass measurements, used in the fit of Ref. [14], were performed with two-photon production, $\gamma\gamma\rightarrow{\eta_{c}}(1S)\rightarrow\mbox{hadrons}$, radiative decays ${{J\mskip-3.0mu/\mskip-2.0mu\psi\mskip 2.0mu}}\rightarrow{\eta_{c}}(1S)\gamma$ and ${\psi{(2S)}}\rightarrow{\eta_{c}}(1S)\gamma$, ${p}{\overline{p}}\rightarrow{\eta_{c}}(1S)\rightarrow\gamma\gamma$, and exclusive $B$ decays, yielding the average value $2983.4\pm 0.5\mathrm{\,Me\kern-1.00006ptV}$. Mass determinations via exclusive $B$ decays, performed at the BaBar and Belle experiments [49, 50, 51], do not provide consistent results. In 2009, the CLEO collaboration observed a significant asymmetry in the line shapes of radiative ${{J\mskip-3.0mu/\mskip-2.0mu\psi\mskip 2.0mu}}\rightarrow\gamma{\eta_{c}}(1S)$ and ${\psi{(2S)}}\rightarrow\gamma{\eta_{c}}(1S)$ transitions [52], which, when ignored, could lead to significant bias in the mass and width measurement via ${J\mskip-3.0mu/\mskip-2.0mu\psi\mskip 2.0mu}$ or $\psi{(2S)}$ radiative decays. Recent BES III results [53, 54] obtained using radiative decays of $\psi{(2S)}$, shifted the world average value by more than two standard deviations. Therefore precise ${\eta_{c}}(1S)$ mass measurements using a different technique are needed. LHCb measured $M_{{\eta_{c}}(1S)}=2982.2\pm 1.5\pm 0.1\mathrm{\,Me\kern-1.00006ptV}$ [15] using ${\eta_{c}}(1S)$ from $b$-hadron decays and reconstructing ${\eta_{c}}(1S)$ via decays to ${p}{\overline{p}}$. A similar situation occurs with the ${\eta_{c}}(1S)$ natural width determination, where recent BES III results obtained using radiative decays of $\psi{(2S)}$ shifted the world average from $29.7\pm 1.0\mathrm{\,Me\kern-1.00006ptV}$ to $31.8\pm 0.8\mathrm{\,Me\kern-1.00006ptV}$.

The properties of the ${\eta_{c}}(2S)$ state are less well studied. Measurements at the CLEO [55], BaBar [56, 57], Belle [58, 51] and BES III [59, 60] experiments, using $\gamma\gamma\rightarrow{\eta_{c}}(2S)\rightarrow\mbox{hadrons}$, double charmonium production in $e^{+}e^{-}$ annihilation, exclusive $B$ decays and radiative transitions of $\psi{(2S)}$, yield the world averages [14] of $3639.4\pm 1.3\mathrm{\,Me\kern-1.00006ptV}$ for the ${\eta_{c}}(2S)$ mass, and $11.3^{+3.2}_{-2.9}\mathrm{\,Me\kern-1.00006ptV}$ for its natural width.

Table 5 presents measurements of the masses of the ${\eta_{c}}$ and $\chi_{c}$ states and of the natural width of the ${\eta_{c}}(1S)$ from the fit of the $\phi\phi$ invariant mass spectrum in Fig. 4. For the determination of the systematic uncertainties, except for the test of the impact of the $f_{0}(980)$ meson, the same variations of the analysis are performed as for the determination of the charmonium yields. In addition, the effect of excluding the ${\eta_{c}}(2S)$ mass region ($2.8-3.7\mathrm{\,Ge\kern-1.00006ptV}$) is studied, and the uncertainties related to the momentum-scale calibration are estimated by varying the calibration parameter by $\pm\,3\times 10^{-4}$ [33]. The resulting total systematic uncertainty is obtained as the quadratic sum of the individual contributions. The uncertainty related to the momentum-scale calibration dominates the mass determination for all ${\eta_{c}}$ and $\chi_{c}$ states. The uncertainty of the $\Gamma_{{\eta_{c}}(1S)}$ measurement is dominated by the background description.

The measured charmonium masses agree with the world averages [14]. The measured ${\eta_{c}}(1S)$ mass is in agreement with the previous LHCb measurement using decays to the ${p}{\overline{p}}$ final states [15] and has a better precision. The precision obtained for the ${\eta_{c}}(1S)$ mass is comparable to the precision of the world average value. The value of the ${\eta_{c}}(1S)$ natural width is consistent with the world average [14].

The charmonium mass differences $M_{{\chi_{{c}1}}}-M_{{\chi_{{c}0}}}$, $M_{{\chi_{{c}2}}}-M_{{\chi_{{c}0}}}$, and $M_{{\eta_{c}}(2S)}-M_{{\eta_{c}}(1S)}$ are obtained (Table 6) as a consistency check and for comparison with theory. For the determination of the systematic uncertainties the same variations of the analysis are performed as for the determination of the charmonium masses and widths. The uncertainty related to the 2D fit dominates the $M_{{\chi_{{c}1}}}-M_{{\chi_{{c}0}}}$ mass difference measurement. The systematic uncertainty of the $M_{{\chi_{{c}2}}}-M_{{\chi_{{c}0}}}$ measurement is dominated by the uncertainty related to potential contributions from other resonances and by the uncertainty on the background model. The uncertainty related to the momentum-scale calibration dominates the $M_{{\eta_{c}}(2S)}-M_{{\eta_{c}}(1S)}$ mass difference measurement. The measured charmonium mass differences agree with the world averages.

Figure 7 shows the $\Gamma_{{\eta_{c}}(1S)},\,M_{{\eta_{c}}(1S)}$ contour plot, obtained from the analysis of $b$-hadron decays into $\eta_{c}$ mesons, where the $\eta_{c}$ candidates are reconstructed via the decay ${\eta_{c}}(1S)\rightarrow\phi\phi$.

The measurements of the ${\eta_{c}}(1S)$ mass and natural width using ${\eta_{c}}(1S)$ meson decays to $\phi\phi$ are consistent with the studies using decays to $p$ $\overline{p}$ [15] and with the world average [14]. The measured ${\eta_{c}}(1S)$ mass is below the result in Ref. [61]. The precision obtained on the ${\eta_{c}}(1S)$ mass is comparable to the precision of the world average.

6 First evidence of the $\boldsymbol{{{B}^{0}_{s}}\rightarrow\phi\phi\phi}$ decay

In order to extract $\phi\phi\phi$ combinations a 3D extended unbinned maximum likelihood fit is used, as described in Sect. 3. Figure 8 shows the invariant mass distribution for $\phi\phi\phi$ combinations.

The fit to the invariant $\phi\phi\phi$ mass spectrum is performed using a sum of two Gaussian functions with a common mean to describe the ${B}^{0}_{s}$ signal, and an exponential function to describe combinatorial background. The ratio of the two Gaussian widths and the fraction of the narrow Gaussian are taken from simulation so that a single free parameter in the $\phi\phi\phi$ invariant mass fit accounts for the detector resolution. A signal of $41\pm 10\pm 5$ ${B}^{0}_{s}$ decays over a low background of about 3 events is obtained. Uncertainties related to the background description in the 3D fit and to the decay model defining the $\phi$ polarization dominate the systematic uncertainty in the ${B}^{0}_{s}$ signal yield determination. The significance of the signal is estimated from the distributions of the difference in the logarithm of the best-fit $\chi^{2}$ with and without including the signal shape in toy simulation samples. This leads to a signal significance of 4.9 standard deviations.

The ${{B}^{0}_{s}}\rightarrow\phi\phi$ decay mode is chosen as a normalization mode for the ${\mathcal{B}}({{B}^{0}_{s}}\rightarrow\phi\phi\phi)$ measurement. The invariant mass spectrum obtained from 2D fits in bins of the $\phi\phi$ invariant mass in the region of the ${B}^{0}_{s}$ mass is shown in Fig. 9.

A sum of two Gaussian functions with a common mean is used to describe the ${B}^{0}_{s}$ signal shape, while an exponential function models the combinatorial background. The ratio of the two Gaussian widths and the fraction of the narrow Gaussian function are taken from simulation. In total $2701\pm 114\pm 84$ ${B}^{0}_{s}$ decays are found. The uncertainties related to the description of the resolution in the 2D fits and the description of the $\phi\phi$ invariant mass resolution dominate the systematic uncertainty in the ${B}^{0}_{s}$ signal yield determination.

The ratio of the ${{B}^{0}_{s}}\rightarrow\phi\phi\phi$ and ${{B}^{0}_{s}}\rightarrow\phi\phi$ branching fractions is obtained from the relative ${{B}^{0}_{s}}\rightarrow\phi\phi\phi$ and ${{B}^{0}_{s}}\rightarrow\phi\phi$ signal yields and their efficiencies as

[TABLE]

In the above expression, the event yields are determined from the fits. The efficiency ratio, $\varepsilon_{{{B}^{0}_{s}}\rightarrow\phi\phi\phi}/\varepsilon_{{{B}^{0}_{s}}\rightarrow\phi\phi}=0.26\pm 0.01$, is obtained from simulation and corrected to account for different ${B}^{0}_{s}$ transverse momentum spectra in data and simulation. The ${{B}^{0}_{s}}\rightarrow\phi\phi\phi$ transition is assumed to proceed via a three-body decay with uniform phase-space density. This assumption is supported below by comparing the $\phi\phi$ invariant mass distribution in data and simulation. The systematic uncertainty is dominated by the uncertainty in polarization of the $\phi$ mesons in the decay ${{B}^{0}_{s}}\rightarrow\phi\phi\phi$, as discussed at the end of this section. Using the branching fraction of the ${{B}^{0}_{s}}\rightarrow\phi\phi$ decay, ${\mathcal{B}}({{B}^{0}_{s}}\rightarrow\phi\phi)=(1.84\pm 0.05\pm 0.07\pm 0.11_{f_{s}/f_{d}}\pm 0.12_{\scalebox{0.6}{norm}})\times 10^{-5}$ [22], the branching fraction for the ${B}^{0}_{s}$ meson decay to three $\phi$ mesons is determined to be

[TABLE]

where the last uncertainty is due to the branching fraction ${\mathcal{B}}({{B}^{0}_{s}}\rightarrow\phi\phi)$.

The ${{B}^{0}_{s}}\rightarrow\phi\phi\phi$ transition can proceed via a two-body decay involving intermediate resonances or via a three-body ${{B}^{0}_{s}}\rightarrow\phi\phi\phi$ decay. In order to search for contributions from possible intermediate resonances, the invariant mass of each $\phi\phi$ combination from all ${{B}^{0}_{s}}\rightarrow\phi\phi\phi$ candidates in the signal region of $\pm 3$ standard deviations around the ${B}^{0}_{s}$ mass is examined, see Fig. 10. The ${B}^{0}_{s}$ candidates are constrained to the known ${B}^{0}_{s}$ mass. Three entries to the histogram are produced by each ${B}^{0}_{s}$ candidate. A phase-space distribution as obtained from simulation is overlaid for comparison. No indication of significant contributions from $\eta_{c}$, ${\chi_{c}}$, $f_{2}(2300)$ or $f_{2}(2340)$ states is seen.

A symmetrized Dalitz plot constructed following the approach described in Ref. [62] shows no evidence for resonant contributions either.

The polarization of the $\phi$ mesons is studied by means of the angle $\theta$ between the direction of flight of a $\phi$ meson in the ${B}^{0}_{s}$ rest frame and the ${B}^{0}_{s}$ direction in the laboratory frame. With the limited sample of ${{B}^{0}_{s}}\rightarrow\phi\phi\phi$ candidates the 3D fit technique to remove contributions from ${{K}^{+}}{{K}^{-}}$ combinations that are not from $\phi$ decays cannot be used for this measurement. Instead, all $\phi$ mesons contributing in the mass range of the ${B}^{0}_{s}$ are used, with an estimated signal purity of $71\%$. Figure 11 compares the $\cos(\theta)$ distribution for the ${{B}^{0}_{s}}\rightarrow\phi\phi\phi$ signal candidates in data with expectations from simulation using different assumptions for the polarization. The purely longitudinal polarization clearly does not describe the data. The difference between the expectations for no polarization and purely transverse polarization is used to estimate the corresponding systematic uncertainty in the ${\mathcal{B}}({{B}^{0}_{s}}\rightarrow\phi\phi\phi)$ measurement. The most probable value for the fraction of transverse polarization, $f_{\textrm{T}}$, is found to be $f_{\textrm{T}}=0.86$. Assuming a uniform prior in the physically allowed range, a Bayesian lower limit of $f_{\textrm{T}}>0.28$ at $95\%$ CL is found.

7 Summary and discussion

Charmonium production in $b$-hadron inclusive decays is studied in $p$ $p$ collisions collected at $\sqrt{s}=7$ and $8\mathrm{\,Te\kern-1.00006ptV}$ corresponding to an integrated luminosity of $3.0\mbox{\,fb}^{-1}$, using charmonium decays to $\phi$-meson pairs. The masses and natural widths of the $\eta_{c}$ and $\chi_{c}$ states are determined. In addition, the first evidence of ${{B}^{0}_{s}}\rightarrow\phi\phi\phi$ decay is obtained.

Ratios of charmonium $C$ production rates,

[TABLE]

are measured to be

[TABLE]

where the first uncertainties are statistical and the second ones are systematic. Using the branching fractions of $\chi_{c}$ decays to $\phi\phi$ from Ref. [14], relative branching fractions of $b$ hadrons decaying inclusively to $\chi_{c}$ states are derived,

[TABLE]

where the third uncertainty is due to the branching fractions ${\mathcal{B}}(\chi_{c}\rightarrow\phi\phi)$. These results are consistent with the ratio of the $\chi_{{c}1}$ and $\chi_{{c}2}$ production rates measured in ${B}^{0}$ and ${{B}^{+}}$ decays [14].

Inclusive production rates of the $\chi_{c}$ states in $b$-hadron decays are derived using branching fractions of the $\chi_{c}$ decays to $\phi\phi$ from Ref. [14], an average of the results from Belle [46] and BaBar [47] ${\mathcal{B}}({\eta_{c}}(1S)\rightarrow\phi\phi)=(3.21\pm 0.72)\times 10^{-3}$, and the ${\eta_{c}}(1S)$ inclusive production rate measured using decays to $p$ $\overline{p}$, ${\mathcal{B}}(b\rightarrow{\eta_{c}}(1S)X)=(4.88\pm 0.97)\times 10^{-3}$ [15]. They are

[TABLE]

where the third uncertainty is due to the uncertainties on the branching fraction of the $b$-hadron decays to the ${\eta_{c}}(1S)$ meson, ${\mathcal{B}}({b}\rightarrow{\eta_{c}}(1S)X)$, and ${\eta_{c}}(1S)$ and $\chi_{c}$ decays to $\phi\phi$. No indirect contribution to the production rate is subtracted. However, since contributions from $\psi{(2S)}$ decays to the $\chi_{c}$ states are limited, the results disfavour dominance of either colour-octet or colour-singlet contributions. The observed relations between the $\chi_{c}$ branching fractions are not consistent with those predicted in Ref. [18]. The branching fraction ${\mathcal{B}}(b\rightarrow{\chi_{{c}0}}X)$ is measured for the first time. The result for $b$-hadron decays into $\chi_{{c}1}$ is the most precise measurement for the mixture of ${B}^{0}$, ${{B}^{+}}$, ${B}^{0}_{s}$, ${B}_{c}^{+}$ and $b$-baryons. The central value of the result for $b$-hadron decays into $\chi_{{c}1}$ is lower than the central values measured by the DELPHI [8] and L3 [9] experiments at LEP. The value obtained is consistent with the branching fraction of $b$-hadron decays into $\chi_{{c}1}$ measured by CLEO [2], Belle [4] and BaBar [5] with the light mixture of ${B}^{0}$ and ${{B}^{+}}$. The branching fraction of $b$-hadron decays into $\chi_{{c}2}$ is measured for the first time with the ${B}^{0}$, ${{B}^{+}}$, ${B}^{0}_{s}$ and $b$-baryons mixture. The result is consistent with the world average corresponding to the ${B}^{0}$, ${{B}^{+}}$ mixture [14] and with individual measurements from CLEO [3], Belle [4], and BaBar [5].

Scaled differential charmonium production cross-sections as a function of $p_{\mathrm{T}}$ are presented for the ${\eta_{c}}(1S)$ and $\chi_{c}$ states in the LHCb acceptance and for \mbox{p_{\mathrm{T}}}>4\mathrm{\,Ge\kern-1.00006ptV}. Next-to-leading-order calculations of the $p_{\mathrm{T}}$ dependence of the $\eta_{c}$ and $\chi_{c}$ production rates in $b$-hadron decays will help to relate the results to conclusions on production mechanisms.

The production rate of the ${\eta_{c}}(2S)$ state in $b$-hadron decays is determined to be

[TABLE]

This is the first measurement for inclusive ${\eta_{c}}(2S)$ production rate in $b$-hadron decays and the first evidence for the decay ${\eta_{c}}(2S)\rightarrow\phi\phi$. The production rate as a function of the assumed natural width is given in Fig. 5. These are the first $\chi_{c}$ and ${\eta_{c}}(2S)$ inclusive production measurements, using charmonium decays to a hadronic final state, in the high-multiplicity environment of a hadron machine. In addition, upper limits at $95\%$ ($90\%$) CL on the production rates of the $X(3872)$, $X(3915)$, and ${\chi_{{c}2}}(2P)$ states in $b$-hadron decays are obtained,

[TABLE]

or

[TABLE]

Masses and natural widths of the $\eta_{c}$ and $\chi_{c}$ states agree with the world averages. The precision of the ${\eta_{c}}(1S)$ mass is comparable to the precision of the world average value. The measured ${\eta_{c}}(1S)$ mass is in agreement with the LHCb measurement using decays to the ${p}{\overline{p}}$ final states [15].

First evidence for the transition ${{B}^{0}_{s}}\rightarrow\phi\phi\phi$ is reported with a significance of 4.9 standard deviations. Its branching fraction is measured to be

[TABLE]

No resonant structure is observed in the $\phi\phi$ invariant mass distribution. In the ${{B}^{0}_{s}}\rightarrow\phi\phi\phi$ decay, transverse polarization is preferred for the $\phi$ mesons, with an estimate of $f_{\textrm{T}}>0.28$ at $95\%$ CL and the most probable value of $f_{\textrm{T}}=0.86$ for the fraction of transverse polarization.

As a by-product of the analysis, the branching fraction ${\mathcal{B}}({{B}^{0}_{s}}\rightarrow\phi\phi)$ is determined to be ${\mathcal{B}}({{B}^{0}_{s}}\rightarrow\phi\phi)=(2.18\pm 0.17\pm 0.11\pm 0.14_{f_{s}}\pm 0.65_{{\mathcal{B}}})\times 10^{-5}$ with a different technique with respect to the previous results [63, 64, 21, 22]. This technique is based on relation of ${B}^{0}_{s}$ production in $p$ $p$ collisions and ${\eta_{c}}(1S)$ inclusive production rate in $b$-hadron decays, and reconstruction of ${B}^{0}_{s}$ and ${\eta_{c}}(1S)$ via decays to $\phi\phi$. The measurement is consistent with the recent LHCb result [22] and the current world average [14], as well as with theoretical calculations [29, 65, 30].

Finally, using the measurements presented and external input, the ratio of the branching fractions for the ${\eta_{c}}(1S)$ decays to $\phi\phi$ and to ${p}{\overline{p}}$ is determined. The measured ${B}^{0}_{s}$ and ${\eta_{c}}(1S)$ yields and efficiency ratio, the branching fraction ${\mathcal{B}}({{B}^{0}_{s}}\rightarrow\phi\phi)=(1.84\pm 0.05\pm 0.07\pm 0.11_{f_{s}/f_{d}}\pm 0.12_{\scalebox{0.6}{norm}})\times 10^{-5}$ [22], the ${J\mskip-3.0mu/\mskip-2.0mu\psi\mskip 2.0mu}$ production rate in $b$-hadron decays ${\mathcal{B}}({b}\rightarrow{{J\mskip-3.0mu/\mskip-2.0mu\psi\mskip 2.0mu}}X)=(1.16\pm 0.10)\%$ [14], the relative production rates of ${\eta_{c}}(1S)$ and ${J\mskip-3.0mu/\mskip-2.0mu\psi\mskip 2.0mu}$ in $b$-hadron decays $\frac{{\mathcal{B}}({b}\rightarrow{\eta_{c}}(1S)X)\times{\mathcal{B}}({\eta_{c}}(1S)\rightarrow{p}{\overline{p}})}{{\mathcal{B}}({b}\rightarrow{{J\mskip-3.0mu/\mskip-2.0mu\psi\mskip 2.0mu}}X)\times{\mathcal{B}}({{J\mskip-3.0mu/\mskip-2.0mu\psi\mskip 2.0mu}}\rightarrow{p}{\overline{p}})}=0.302\pm 0.042$ [15], the branching fraction ${\mathcal{B}}({{J\mskip-3.0mu/\mskip-2.0mu\psi\mskip 2.0mu}}\rightarrow{p}{\overline{p}})=(2.120\pm 0.029)\times 10^{-3}$ [14], the ratio of fragmentation fractions $f_{s}/f_{d}=0.259\pm 0.015$ [66], and the ${\mathchar 28931\relax}^{0}_{b}$ fragmentation fraction $f_{{{\mathchar 28931\relax}^{0}_{b}}}$ momentum dependence from Ref. [67] are used. The ratio of the branching fractions for the ${\eta_{c}}(1S)$ decays to $\phi\phi$ and to ${p}{\overline{p}}$ is determined as

[TABLE]

where the third uncertainty is related to $f_{s}/f_{d}$, the fourth uncertainty is related to $f_{{{\mathchar 28931\relax}^{0}_{b}}}$, and the fifth uncertainty is related to uncertainties of the production rates and decay branching fractions involved. This value is larger than the value computed from the world average branching fractions given in Ref. [14].

Acknowledgements

We would like to thank Emi Kou for motivating the studies of charmonium production in LHCb using hadronic final states and the useful discussions regarding charmonium production mechanisms. 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 (The Netherlands); MNiSW and NCN (Poland); MEN/IFA (Romania); MinES and FASO (Russia); MinECo (Spain); SNSF and SER (Switzerland); NASU (Ukraine); STFC (United Kingdom); NSF (USA). We acknowledge the computing resources that are provided by CERN, IN2P3 (France), KIT and DESY (Germany), INFN (Italy), SURF (The 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łodowska-Curie Actions and ERC (European Union), Conseil Général de Haute-Savoie, Labex ENIGMASS and OCEVU, Région Auvergne (France), RFBR and Yandex LLC (Russia), GVA, XuntaGal and GENCAT (Spain), Herchel Smith Fund, The Royal Society, Royal Commission for the Exhibition of 1851 and the Leverhulme Trust (United Kingdom).