1 Introduction

The constituent quark model [1,2,3] has been very successful describing the quark composition and allowed quantum numbers of observed hadron states. Potential models exploiting heavy-quark symmetry [4,5,6] are used to calculate properties of mesons containing one heavy and one light quark, including those with \(\bar{b}s\) quark content, termed collectively the \(B_s^{**0}\) mesons. Yet it is still difficult to predict precise masses and widths of such states. Experimental measurements are therefore essential. This paper presents an observation of an excess in the \({{B} ^+} {{K} ^-} \) invariant-mass spectrum,Footnote 1 which we interpret as originating from \(B_s^{**0}\) mesons. Two states, \(B_{s1}^0\) and \(B_{s2}^{*0}\), with masses higher than the ground state \({B} ^0_{s} \) and \(B_s^{*0}\) mesons have been previously observed [7,8,9,10,11]. Based on their masses and narrow widths, they are interpreted as two of the predicted states with orbital angular momentum \(L=1\). We report the ratio of the production cross-section times branching fraction for the newly observed states relative to that of the \(B_{s2}^{*0}\) meson.

Predictions using a variety of methods have been made for masses and widths of states belonging to orbital or radial excitations of \(\bar{b}s\) states [12,13,14,15,16,17,18]. The spectroscopic notation, \(n^{2S+1}L_J\), is commonly used to refer to these states, where n gives the radial quantum number, S the sum of the quark spins, L the orbital angular momentum of the quarks, and J the total angular momentum. The masses of the first radially excited states, \(2^1S_0\) and \(2^3S_1\), are predicted in the range between 5920\(\,\text {Me}\)V and 6020\(\,\text {Me}\)V.Footnote 2 Different calculations give significantly different predicted widths for these states. In some cases they are close to 100\(\,\text {Me}\)V or more [12,13,14] and in others can be less than 50\(\,\text {Me}\)V  [15,16,17]. Broad states are difficult to identify experimentally, however, because of large non-resonant continuum contributions. The first D-wave states (\(L=2\)) are predicted to have masses ranging roughly from 6050 to 6200\(\,\text {Me}\)V. The \(1^3D_3\) state is almost always predicted to be relatively narrow, with a width less than 50\(\,\text {Me}\)V. The \(J=2\) states \(1^1D_2\) and \(1^3D_2\) can mix, and depending on the mixing angle they may produce one of the states as narrow as 20\(\,\text {Me}\)V. With multiple states potentially decaying to both the \({{B} ^+} {{K} ^-} \) final state directly and through a \(B^{*+}\) intermediate state, with a missing photon from the \(B^{*+} \rightarrow B^+ \gamma \) decay, the observation of overlapping peaks in this mass region seems likely.

2 Detector and selection

The LHCb detector [19, 20] 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 momentum to 1.0% at 200\(\,\text {Ge}\)V. The minimum distance of a track to a primary pp interaction vertex (PV), the impact parameter (IP), is measured with a resolution of \((15+29/p_{\mathrm {T}})\,\upmu \text {m} \), where \(p_{\mathrm {T}}\) is the component of the momentum transverse to the beam, in \(\,\text {Ge}\)V. 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. At the hardware trigger stage, events are required to have a muon with high \(p_{\mathrm {T}}\) or a hadron with high transverse energy in the calorimeter. The software trigger requires a two-, three- or four-track vertex with a significant displacement from any PV. At least one charged particle must have a \(p_{\mathrm {T}} > 1.6\,\text {Ge}\text {V} \) and be inconsistent with originating from a PV. Consistency with a PV is defined based on the \(\chi ^2\) difference between vertex fits including and excluding the particle under consideration.

We use data samples collected from 2011 to 2018, at centre-of-mass energies of 7, 8, and 13 TeV, corresponding to an integrated luminosity of 9\(\,\text {fb} ^{-1}\). We simulate both \(B_{s2}^{*0}\) decays and decays of a \(B_s^{**0}\) meson with a mass of 300\(\,\text {Me}\)V above the \({B} ^+\) \({K} ^-\) mass threshold. In the simulation, pp collisions are generated using Pythia  [21, 22] with a specific LHCb configuration [23]. Decays of hadronic particles are described by EvtGen  [24]. The interaction of the generated particles with the detector and its response are implemented using the Geant4 toolkit [25, 26] as described in Ref. [27].

The selection of the \({B} ^+\) \({K} ^-\) candidates is performed in two steps. First, we select \({B} ^+\) candidates using the decays \({{B} ^+} \rightarrow {{J /\psi }} {{K} ^+} \) and \({{B} ^+} \rightarrow {{\overline{D}} {}^0} {{\pi } ^+} \), with \({{J /\psi }} \rightarrow {\mu ^+} {\mu ^-} \) and \({{\overline{D}} {}^0} \rightarrow {{K} ^+} {{\pi } ^-} \). Each final-state particle is formed from a high-quality track that is inconsistent with being produced at any PV in the event. Each kaon, pion and muon is also required to have particle identification information consistent with its hypothesis. We do not require the candidate to have fired the trigger in the event. The selection of \({B} ^+\) candidates has high purity, with a background contribution in the \({B} ^+\) mass region of less than 10%.

Subsequently, we form the \(B_s^{**0}\) candidates by combining the \({B} ^+\) candidate with a \(K^-\) candidate consistent with being produced at the PV, referred to as the prompt kaon. In the same way we also select a separate sample of same-sign \({B} ^+\) \({K} ^+\) candidates further used for constraining combinatorial background. Strong particle-identification requirements are applied to the prompt kaons to reduce the large pion background. As observed for the known \(B_{s2}^{*0}\) resonance, the strongest discriminant between resonant signals and combinatorial background is the kaon transverse momentum. We therefore study the \(B_s^{**0}\) candidates in bins of the prompt kaon \(p_{\mathrm {T}}\): \(0.5< p_{\mathrm {T}} < {1}~\hbox {GeV}\), \(1< p_{\mathrm {T}} < {2}~\hbox {GeV}\), and \(p_{\mathrm {T}} > {2}~\hbox {GeV}\). The \(B_s^{**0}\) candidate mass, \(m_{{{B} ^+} {{K} ^-}}\), is reconstructed constraining the masses \({J /\psi }\) or \({\overline{D}} {}^0\) meson and the \({B} ^+\) meson to their known values [28] and requiring the \({B} ^+\) candidate to have been produced at the PV. We search for new states in a spectrum of the mass difference, \(\Delta m = m_{{{B} ^+} {{K} ^-}} - m_{{{B} ^+}} - m_{{{K} ^-}}\), where the latter two masses refer to the known \({{B} ^+} \) and \({K} ^-\) masses [28]. The \(\Delta m\) spectrum, measured in the three kaon \(p_{\mathrm {T}}\) regions, is shown in Fig. 1. A clear deviation from a smooth continuum shape is observed at approximately 300\(\,\text {Me}\)V above the mass threshold, and is consistently present in both \({B} ^+\) selections.

3 Fit description and results

We determine the significance of the excess using an unbinned maximum-likelihood fit to the \({B} ^+\) \({K} ^-\) mass difference spectrum in the region from 150\(\,\text {Me}\)V to 800 MeV above the kinematic threshold. The fit is performed simultaneously in three bins of the prompt kaon \(p_{\mathrm {T}}\).

We form the background distribution from two components. The first is the combinatorial background, whose shape and yield are fixed based on a fit to the same-sign \({B} ^+\) \({K} ^+\) mass difference distribution using a threshold function of the form \(f(\Delta m) = \Delta m^{a}\exp (-b\Delta m)\), with a and b as free parameters. The additional continuum background component that is not present in the same-sign sample, and which arises from associated production (AP) of the non-resonant \({B} ^+\) \({K} ^-\) pairs, is free to vary in the fit; our model for this component is a second-degree polynomial. Given the similarity of the combinatorial and AP background shapes, statistical uncertainties in the former are expressed through the latter fit component.

The \(B_s^{**0}\) signal is described by one or two relativistic Breit–Wigner distributions, convolved with resolution functions determined using simulation. For new states, the decay products’ angular momentum (\(\ell \)) is unknown, making it difficult to assign an appropriate Blatt–Weisskopf barrier factor [29]. The \(B_s^{**0}\) state most consistently predicted to be narrow decays with a \(B^{(*)+}{{K} ^-} \) angular momentum of \(\ell =3\), so we use this as default, but consider other angular momenta in alternative fits. The barrier radius parameter is fixed to \({4}~\hbox {GeV}^{-1}\) [30, 31]. When fitting with two peaks, no interference between the two new states is considered. Interference is possible if the two states are both \(J=2\) mixed D-wave states, but in that case one of the states should be much wider, so this scenario is unlikely to produce two narrow peaks.

We use two resolution models: one where the decay goes directly to a \({B} ^+\) \({K} ^-\) final state and the other where it proceeds through an intermediate \(B^{*+}\) meson which decays to \({{B} ^+} \gamma \). We model the resolution in the direct \({B} ^+\) \({K} ^-\) decay as the sum of a Gaussian distribution and a Crystal Ball function [32] with a shared mean, and in the \(B^{*+}{{K} ^-} \) channel as the sum of two Gaussian distributions with different means to account for the smearing and shift in the peak position by approximately 45 MeV due to the unreconstructed photon. The detector resolution is small compared to the width of the prospective states; the effect of the missing photon is significant. The peak positions and widths are shared in each bin of the kaon \(p_{\mathrm {T}}\), but the yield is allowed to vary independently. The fit results using the polynomial AP shape are shown in Fig. 1. A total signal yield of \({18\,900 \pm 2200}\) is observed in the two peaks across all the \(p_{\mathrm {T}}\) bins when fit with the default model and \({B} ^+\) \({K} ^-\) resolution.

Fig. 1
figure 1

The \({B} ^+\) \({K} ^-\) mass difference distributions in data, overlaid with the fit models: (left column) one-peak hypothesis and (right column) two-peak hypothesis. In each column, the rows are from top to bottom for candidates with the prompt kaon \(p_{\mathrm {T}}\): \(0.5< p_{\mathrm {T}} < {1}~\hbox {GeV}\), \(1< p_{\mathrm {T}} < {2}~\hbox {GeV}\), and \(p_{\mathrm {T}} > {2}~\hbox {GeV}\). The legend in the top row applies for each column. The associated production (AP) background is described by a second-degree polynomial in each fit. The combinatorial background shape is fixed from a fit to the \({B} ^+\) \({K} ^+\) mass difference distributions

Full size image

We determine the local statistical significance of the observation using the asymptotic formula for the profile likelihood ratio test statistic [33]. We compare first the one-peak hypothesis to the null hypothesis, and then the two-peak to the one-peak hypothesis. The null fit model corresponds to only a polynomial description of the AP in addition to the shape of the same-sign data. The minimum of the test statistics across each fit with different AP descriptions (discussed later) gives local significances of more than \(20\sigma \) for the one-peak fit with respect to the null hypothesis and \(7.7\sigma \) for the two-peak description with respect to the one-peak hypothesis. As the significance is large, we have not made an exact quantification of the penalty due to the look-elsewhere effect.

Assuming the two-peak hypothesis, the masses and widths of the two states and the fraction of the yield in the lower mass state, \(f_1\), are determined to be

$$\begin{aligned} m_1&= 6063.5 \pm 1.2 \text { (stat)} \pm 0.8\text { (syst)}\,\text {Me}\text {V}, \\ \Gamma _1&= 26 \pm 4 \text { (stat)} \pm 4\text { (syst)}\,\text {Me}\text {V}, \\ m_2&= 6114 \pm 3 \text { (stat)} \pm 5\text { (syst)}\,\text {Me}\text {V}, \\ \Gamma _2&= 66 \pm 18 \text { (stat)} \pm 21\text { (syst)}\,\text {Me}\text {V}, \\ f_1&= 0.47 \pm 0.11 \text { (stat)} \pm 0.15\text { (syst)}, \end{aligned}$$

if the decay is directly to the \({B} ^+\) \({K} ^-\) final state. If it instead proceeds through \(B^{*+}{{K} ^-} \), they are

$$\begin{aligned} m_1&= 6108.8 \pm 1.1\text { (stat)} \pm 0.7\text { (syst)}\,\text {Me}\text {V}, \\ \Gamma _1&= 22 \pm 5\text { (stat)} \pm 4\text { (syst)}\,\text {Me}\text {V}, \\ m_2&= 6158 \pm 4\text { (stat)} \pm 5\text { (syst)}\,\text {Me}\text {V}, \\ \Gamma _2&= 72 \pm 18\text { (stat)} \pm 25\text { (syst)}\,\text {Me}\text {V}, \\ f_1&= 0.42 \pm 0.11\text { (stat)} \pm 0.16\text { (syst)}. \end{aligned}$$

The dominant contribution to the systematic uncertainty arises from the considered variations in the fit functions, and is obtained by repeating the fits with different conditions and calculating the differences with the default fit. We refit the data sample using alternative functional forms for the AP background: an exponential and a separate threshold function. Each AP description produces consistent results for the significance, with variation in the peak parameters. The dependence on the signal model is estimated by alternatively using relativistic Breit-Wigner models with different angular-momentum hypotheses of \(\ell =2\) or \(\ell =1\) and taking the largest absolute difference. We also vary the \(p_{\mathrm {T}}\) binning and mass range for the fit, and we include the effect on one peak’s parameters of changing the other peak from the \({B} ^+\) \({K} ^-\) model to the \(B^{*+}{{K} ^-} \) model or vice versa. A breakdown of the systematic uncertainties is given in Table 1.

Table 1 Sources of systematic uncertainty on the peak parameters assuming a two-peak structure decaying either directly to \({{B} ^+} {{K} ^-} \) or through \(B^{*+}{{K} ^-} \). The parameter differences with respect to the nominal fit result are given in \(\,\text {Me}\)V
Full size table

Since the mass difference between the two peaks is found to be close to the \(B^{*+}\)–\({B} ^+\) mass difference of approximately 45 MeV, we consider also the possibility that the structure is produced by a single resonance that decays in both the \({B} ^+\) \({K} ^-\) and \(B^{*+}{{K} ^-} \) channels. In an alternative fit, both resolution models are simultaneously combined with a resonance lineshape described by a single mass and width. This hypothesis is disfavored by more than \(2\sigma \) with respect to the two-state hypothesis and cannot be completely ruled out. Qualitatively, this model does not describe the data as well because the observed width is much different in the two peaks, with the lower peak being narrower.

4 Production ratio

In addition to the masses and widths of the new mesons, we also determine the production ratio relative to the \(B_{s2}^{*0}\) meson. The \({B} ^+\) \({K} ^-\) candidates in both mass ranges are selected as previously described, but we use only candidates selected from the \({{B} ^+} \rightarrow {{J /\psi }} {{K} ^+} \) decay and require the candidate to have triggered the event. We define the production ratio as the product of the cross-sections times branching fractions of the new states divided by the corresponding product for \(B_{s2}^{*0}\),

$$\begin{aligned} R = \frac{ \sum \sigma ({B_s^{**0}}) \times {\mathcal {B}} ( {B_s^{**0}} \rightarrow {B} ^{(*)+}{{K} ^-} ) }{ \sigma ({B_{s2}^{*0}})\times {\mathcal {B}} ({B_{s2}^{*0}} \rightarrow {{B} ^+} {{K} ^-} ) }. \end{aligned}$$
(1)

We sum the measured yield for both observed states, and consider the yield difference of the \({{B} ^+} {{K} ^-} \) and \({B^{*+}} {{K} ^-} \) models as a systematic uncertainty. The ratio is taken with respect to the yield for only the \({B} ^+\) \({K} ^-\) decay of the \(B_{s2}^{*0}\) meson. The production cross-section is restricted to the fiducial region of \(B_s^{**0}\) transverse momentum and rapidity, \(0< p_{\mathrm {T}} < {20}~\hbox {GeV}\) and \(2.1< y < 4.2\).

The relative efficiency for a state at the observed mass with respect to the \(B_{s2}^{*0}\) meson is determined using simulation, and used to correct the observed yields. The largest systematic uncertainty on the relative efficiency comes from the unknown transverse momentum spectra of the new states. Additional efficiency uncertainties, for example from the description of particle identification in simulation, give smaller contributions. The \(B_{s2}^{*0}\) yield is determined from a separate unbinned maximum-likelihood fit in the range \(50< \Delta m < {85}~\hbox {MeV}\), using the same prompt kaon \(p_{\mathrm {T}}\) binning as in the main fit. The systematic uncertainties from signal and background fit-model variations are handled as in the fit used for the observation. The fit is additionally separated into different data-taking periods since the efficiency and production ratio may depend on the collision energy and detector conditions. The results, however, are statistically compatible between the different periods. We therefore quote the combined production ratio

$$\begin{aligned} R = 0.87 \pm 0.15 \text { (stat)} \pm 0.19 \text { (syst)}. \end{aligned}$$

A breakdown of the systematic uncertainties is given in Table 2.

Table 2 Summary of the systematic uncertainties on the production ratio. The effect of the unknown \(p_{\mathrm {T}}\) distribution is separated from other uncertainties affecting the relative efficiency
Full size table

5 Conclusions

An excess is observed in the \({B} ^+\) \({K} ^-\) mass spectrum at a mass approximately 300\(\,\text {Me}\)V above the \({B} ^+\) \({K} ^-\) threshold, which we interpret as resulting from the decays of multiple excited \(B_s^{**0}\) mesons. The structure is not well described by a single resonance; the significance for a two-peak structure relative to a one-peak hypothesis is greater than \(7\sigma \). A single resonance which can decay through both the \({B} ^+\) \({K} ^-\) and \(B^{*+}\) \({K} ^-\) channels is disfavoured but cannot be excluded. Future investigation with a larger data sample and with the \(B^{*+}\) \({K} ^-\) final state fully reconstructed could potentially better determine the exact structure in this mass range. Under the two-peak hypothesis, we have measured the masses and widths of the two states. Based on predictions for the masses of excited \(B_s^{*0}\) states, the new observation is likely to result from \(L=2\) orbitally excited mesons. We have additionally determined the production cross-section times branching fraction of the new excess relative to the previously observed \(L=1\) \(B_{s2}^{*0}\), which shows that it is produced at a comparable rate.