1 Introduction

Light GeV-scale Higgs-like scalars \(\varphi \) occur in several well-motivated extensions of the Standard Model (SM), in which a SM gauge group singlet scalar field \(\varphi \) couples to the Higgs via a cubic or quartic interaction. Examples include a mediator between a dark sector and the SM [1, 2], a light inflaton [3,4,5], the pseudo-Goldstone bosons associated with the breaking of scale-invariance in a classically scale-invariant model [6], or light scalars to explain different low-energy anomalies, see e.g. [7,8,9]. The light scalar could be the scalar breaking \(B-L\) symmetry [10, 11], which has been studied in e.g. [12].

The phenomenology of GeV-scale Higgs-like scalars has been recently studied in [13,14,15] (see also the review [16] and [17] for recent stellar limits). Below the scale of the electroweak symmetry breaking, the interaction of scalars with the SM particles may be generically described by the two independent interactions: the mixing with the Higgs boson h (parameterised by the mixing angle \(\theta \)), and the trilinear coupling \(h \varphi \varphi \) which gives rise to an invisible Higgs decay. The mixing coupling makes the scalar unstable: it may decay into a pair of leptons or into hadrons.

The most stringent constraints on Higgs-like scalars with masses \(0.3\,\textrm{GeV}\lesssim m_\varphi <m_B-m_K\) (see e.g. [18]) come from the LHCb displaced vertex search for \(B\rightarrow K \varphi (\rightarrow \mu \bar{\mu })\) [19, 20] which constrains the scalar mixing down to \(\theta \simeq 10^{-4}\) depending on the scalar mass. The corresponding search for long-lived particles decaying to a pair of light leptons, \(e\bar{e}, \mu \bar{\mu }\) or light mesons \(\pi \pi , KK\) [21, 22] are currently less sensitive. However, these constraints are subject to the assumption that the scalar dominantly decays visibly into SM particles which may not hold in general (if additional couplings apart from the mixing exist). Thus it is crucial to also consider searches for invisible decays of the scalar.

Table 1 Upper bounds on the branching ratios \(B\rightarrow K+\textrm{inv}\) and the baseline (improved) expected sensitivities of Belle II to the signal strength relative to the SM assuming the kinematic distribution of a 3-body decay \(B\rightarrow K\nu \bar{\nu }\). The SM prediction has been obtained using the expressions in [27] with the form factors [28] using the Bharucha–Straub–Zwicky (BSZ) parametrization [29] and the recent lattice result for \(B\rightarrow K\) [30]. In addition, we add a 20% correction for \(B\rightarrow K^* \nu \bar{\nu }\) due to finite width effects of \(K^*\) [31]. They agree within errors with [32], which predicts a slightly enhanced branching ratio for \(B^+\rightarrow K^+\nu \bar{\nu }\) and a reduced one for \(B^0\rightarrow K^{*0}\nu \bar{\nu }\). Compared to [33,34,35,36,37], the branching ratios for \(B\rightarrow K \nu \bar{\nu }\) are 5 % larger and the ones for \(B\rightarrow K^* \nu \bar{\nu }\) 10% smaller due to the difference in form factors

The B factories BaBar, Belle and Belle II searched for \(B\rightarrow K^{(*)} +\textrm{inv}\) and reported upper limits on the different decay channels which are listed in Table 1. Belle II is expected to measure all four decay channels including the polarisation fraction \(F_L\) of the \(K^*\) [38], which may be used to search for additional invisible final states like heavy neutral leptons (HNLs) [27, 39] or light invisibly-decaying scalars [35, 40, 41], which may constitute dark matter. In fact, the Belle II collaboration already published their first analysis [42]. A simple weighted average indicates the branching ratio \(\textrm{Br} (B^+\rightarrow K^+ +\textrm{inv}) = (1.1\pm 0.4)\times 10^{-5}\), which is \(1.4\sigma \) in excess over the SM prediction.

The complementarity of displaced and invisible searches has also recently been highlighted for Belle II. Ref. [43] studied a Higgs-like scalar decaying to dark sector particles. It stressed the benefit of searching for displaced pairs of leptons or light mesons in the Belle II detector. They conclude that Belle II is able to probe mixing angles down to \(10^{-5}\). Ref. [44] studied axion-like particles. While displaced vertex searches are currently more sensitive for masses above the muon threshold, invisible decay searches are more sensitive for lighter masses.

In this paper, we study the bounds on the light scalar mixing with the SM Higgs under the assumption of a sizeable invisible width of the scalar \(\varGamma _\textrm{inv}\). We show how the limits from \(B\rightarrow K\mu \bar{\mu }\) become weaker if the scalar invisible width is increased. At the same time bounds from \(B\rightarrow K +\textrm{inv}\) become relevant and for sufficiently large values of \(\varGamma _\textrm{inv}\) they will dominate the bound on the scalar/Higgs mixing. We illustrate this complementarity by using current LHCb bounds on \(B\rightarrow K\mu \bar{\mu }\) [19, 20] and the BaBar bound on \(\textrm{Br}(B^+\rightarrow K^+ +\textrm{inv})\) [24]. We will give also two simple examples for light new physics, where the required values for \(\varGamma _\textrm{inv}\) can be achieved by letting the scalar decay into HNLs.

The paper is structured as follows. In Sect. 2 we introduce the couplings of the scalar and the relevant decay modes. In Sect. 3 we provide an overview of searches for scalars at B factories. The complementarity of \(B\rightarrow K \mu \bar{\mu }\) and \(B\rightarrow K +\textrm{inv}\) is presented as the main result in Sect. 4. In Sect. 5 we provide a few examples for models with a sizeable invisible decay width, before concluding in Sect. 6. Technical details are collected in the appendix.

2 Phenomenology of GeV-scale Higgs-like scalars in a nutshell

We consider a light real scalar fieldFootnote 1\(\varphi \), which mixes with the Higgs with mixing angle \(\theta \). In the minimal scenario, \(\theta \) and the scalar mass \(m_{\varphi }\) control every observable such as the scalar lifetime \(\tau _{S}\propto f(m_{\varphi })\theta ^{-2}\) and the partial branching ratios.

Because of the mixing, the structure of the scalar interaction with SM particles is similar to those for the Higgs, but the mixing angle suppresses the couplings. This way, \(\varphi \) couples to all SM fermions at tree level proportional to their respective masses as \((m_f \sin \theta /v)\). These tree-level interactions generate effective couplings to other SM particles such as gluons, photons, and the bound states such as nucleons [15]. They differ from those of the Higgs boson not only by \(\theta \) but also by mass which determines the scale associated with the couplings.

Thus, the scalar decays to kinematically-accessible SM lepton pairs with partial decay widthFootnote 2

$$\begin{aligned} \varGamma (\varphi \rightarrow f\bar{f})= \frac{\sqrt{2}G_F m_f^2 m_\varphi \sin ^2\theta }{8\pi } \left( 1-\frac{4m_f^2}{m_\varphi ^2}\right) ^{3/2}. \end{aligned}$$
(1)

Decays into hadrons are more complicated. Their description depends on the scalar mass which sets a characteristic energy scale of the process. For \(m_{\varphi }\gg \varLambda _{\textrm{QCD}} ={\mathcal {O}}(1\,\text {GeV})\), they may be described inclusively by the decay of the scalar into a quark–antiquark pair; the corresponding decay width is given by Eq. (1) with the additional colour factor \(N_c=3\). For lower masses, \(m_{\varphi }\simeq \varLambda _{\textrm{QCD}}\), the perturbative QCD description breaks down, and hadronic decays must be treated exclusively, i.e., into different mesons.

The lightest possible hadronic decay is into a pair of pions, \(\varphi \rightarrow \pi \pi \). It may be described using chiral perturbation theory, which, however, quickly becomes unreliable just above the decay threshold due to strong interactions of the pions [45]. Alternatively, the calculation of the decay form-factor may be performed using the method of dispersion relations, see [46] for an overview. It was realised (see, e.g., [47]) that the approach suffers from theoretical uncertainties that may significantly change the results. Recently, [48] has calculated the decay width using experimental data for the gravitational pion form-factors with an uncertainty of about a factor \({\mathcal {O}}(2)\).

Fig. 1
figure 1

Decays of a GeV-scale scalar \(\varphi \) in the minimal model described by the two parameters – the scalar mass \(m_{\varphi }\) and the mixing angle with the SM Higgs boson \(\theta \). The left panel shows the decay width of the scalar at \(\sin ^2\theta = 10^{-4}\), while the right panel the branching ratio of its decays into pairs of muons or hadrons. The blue lines indicate the result obtained in [46], with included NLO corrections for decays into quarks and gluons. The red lines show the calculations of [48]. The difference between them indicates the uncertainty in the description of hadronic decays of the scalar (see text for details)

Other hadronic decays of \(\varphi \), e.g., into a pair of kaons and multihadronic states, also suffer from similar problems. In particular, there is no clear way to describe the transition between the exclusive and inclusive approaches. Refs. [46, 48] match the two regimes at different scalar masses – \(2\,\text {GeV}\) and \(1.3\,\text {GeV}\), respectively. Motivated by this issue, in Fig. 1 we show the lifetime and partial branching ratios of \(\varphi \) following the descriptions from [46] and [48], interpreting the difference between them as the uncertainty in the decay width. Depending on the scalar mass, the difference between the widths may be as large as an order of magnitude. We will comment on the impact of this uncertainty on searches for scalars considered in this work in Sect. 4.

Let us now discuss the scalar production channels. Depending on the facility, the main channels are [15] proton bremsstrahlung or decays of various particles: 2- and 3-body decays of kaons and B mesons \(B\rightarrow \varphi + \text {other}\), or a 2-body decay of the Higgs boson \(h\rightarrow \varphi \varphi \). The production channel of main interest in this work is the decay of B mesons into a scalar and a strange meson:

$$\begin{aligned}&B\rightarrow \varphi +X_{s},\nonumber \\&\quad X_{s} = K,K^{*},K_{0}^{*},K_{1}(1270), K_{1}(1400),K_{2}^{*},\dots \end{aligned}$$
(2)

The decay vertex originates from a flavour-changing neutral current coming from electroweak 1-loop contributions [49, 50] (see also Appendix A for details).

In the limit of small scalar masses \(m_{\varphi }\ll m_{B}\) and \(\sin ^{2}\theta \ll 1\), the total exclusive branching ratio is [15] (see also Appendix A)

$$\begin{aligned} \sum _{X_{s}}\text {Br}(B\rightarrow X_{s}+\varphi ) \approx 3.3\, \sin ^{2}\theta \end{aligned}$$
(3)

Decays into K consist of only around \(10\%\) of this number:

$$\begin{aligned} \text {Br}(B\rightarrow K+\varphi )\approx 0.4\,\sin ^{2}\theta , \quad m_{\varphi }\ll m_{B} \end{aligned}$$
(4)

Despite this small branching ratio, this channel is attractive from the experimental point of view. The main reason is that, unlike the heavier resonances \((K^*, K^*_0, K_1,\ldots )\), K is stable on the experiment scales, so the kinematics reconstruction is simple: it only requires to reconstruct the kaon itself. The other resonances are short-lived and decay into a kaon plus other particles such as pions and photons, which require additional signal selection and reconstruction. Because of this, searches for new physics via the \(B\rightarrow K\mu \bar{\mu }\) channel typically give the strongest constraint.

Fig. 2
figure 2

Left panel: constraints on the scalar parameter space coming from searches for \(B\rightarrow K \mu \bar{\mu }\) at LHCb for \(\varGamma _\textrm{inv} = 0\) (see text for details). Right panel: the scalar lifetime at \(\sin \theta = 10^{-4}\) as a function of its mass for various values of the decay width \(\varGamma _\textrm{inv}\), which, being summed with the \(\theta \)-controlled width gives the total decay width of the scalar. The dashed grey line denotes the lifetime \(\tau _{\phi } = 10\,\text {ps}\), above which the scalar lifetime is too large to decay inside the detector frequently, such that the decay probability scales as \(\tau _{\varphi }^{-1}\) and the event rate becomes independent of the total decay width

3 Overview of searches for scalars at B factories

Two types of searches at B factories – BaBar, Belle, Belle II, and LHCb – may be applied to dark scalars: \(B\rightarrow K+\text {visible}\), or \(B\rightarrow K+\text {invisible}\).

The first type corresponds to the production \(B\rightarrow K+\varphi \), with \(\varphi \) subsequently decaying into visible particles within the detector. In general, such decays may be \(\varphi \rightarrow \mu \bar{\mu }, \varphi \rightarrow \pi \pi , \varphi \rightarrow KK\), as well as decays into jets – quarks and gluons. The number of events for this signature (without taking into account further selection and reconstruction) is proportional to

$$\begin{aligned} \textrm{BR}(B^+\rightarrow K^+ \varphi ) \,\textrm{BR}(\varphi \rightarrow \text {visible}) \, \left[ 1- P(r_\textrm{det}|\beta \gamma )\right] , \end{aligned}$$
(5)

where \(P(r_\textrm{det}|\beta \gamma )=\exp [-r_\textrm{det}/\beta \gamma c\tau _{\varphi }]\) denotes the probability to decay outside of the detector of size \(r_\textrm{det}\)Footnote 3, and \(\textrm{BR}(\varphi \rightarrow \text {visible})\) is the branching ratio of decays of \(\varphi \) into visible particles. The scalar \(\varphi \) is produced on-shell and thus the process can be considered as a series of 2-body decays with

$$\begin{aligned}&\varGamma (B\rightarrow K^{(*)} \varphi (\rightarrow f\bar{f})) \nonumber \\&\quad = \int _0^\infty \frac{dq^2}{2\pi } \frac{\left[ \varGamma (B\rightarrow K^{(*)}\varphi )2m_\varphi \varGamma (\varphi \rightarrow f\bar{f}) \right] _{m_\varphi ^2 \rightarrow q^2}}{(q^2-m_\varphi ^2)^2 + m_\varphi ^2\varGamma ^2_\varphi }\nonumber \\&\quad {\mathop {\longrightarrow }\limits ^{\varGamma _\varphi \ll m_\varphi }} \varGamma (B\rightarrow K^{(*)} \varphi ) \textrm{BR}(\varphi \rightarrow f\bar{f}) \end{aligned}$$
(6)

where \(q^2\) denotes the 4-momentum squared of the scalar and the narrow width approximation has been employed in the last line, which is a good approximation for the relevant parameter space. As a result, the invariant mass distribution of the visible particles (if they are fully reconstructed) would have a peak at \(m_\textrm{inv} = m_{\varphi }\) with a width due to the finite resolution of the 4-momenta reconstruction, which may be used to reduce the background efficiently.

From Fig. 1, we conclude that the most probable decay of a GeV-scale scalar is into hadrons, in particular into a pair of pions, kaons, or jets such as two gluons. The decay into muons is significantly suppressed. However, the latter may give the cleanest decay due to better reconstruction capabilities for muons and lower backgrounds.

For the minimal scalar model with only two parameters \(m_{\varphi }\) and \(\theta \), the most stringent constraints on a GeV-scale Higgs-like scalar come from LHCb [19, 20]. Namely, Ref. [20] searched for displaced vertices in the decays \(B^+\rightarrow K^+ \varphi (\rightarrow \mu \bar{\mu })\) [20], while the work [19] considered \(B^0 \rightarrow K^{*0} \varphi (\rightarrow \mu \bar{\mu })\). Belle II [22] also searched for \(B^+\rightarrow K^+ \varphi (\rightarrow \mu \bar{\mu })\). The strongest constraint is placed by [20] which we thus use it in the following analysis. It placed constraints on the scalar mixing angle as a function of the scalar mass \(m_\varphi \) and the lifetime \(\tau _\varphi \) down to \(\sin \theta \lesssim 10^{-4}\).

We extracted the constraint on BR(\(B^+\rightarrow K^+ \varphi \)) \(\times \) BR(\(\varphi \rightarrow \mu \bar{\mu }\)) as a function of the scalar mass \(m_\varphi \) and lifetime \(\tau _\varphi \) from Fig. 4 in [20] which is model-independent. In Fig. 2 (left) we use this constraint to set a limit on the scalar mixing, where we assume the scalar decay width as in [46] but including NLO corrections to decay widths into quarks and gluons. The gaps at several masses are caused by the exclusions in the searched mass range due to the contribution of SM resonances \(K^{0}_{S}, J/\psi ,\psi (2S)\), and \(\psi (3770)\). Interestingly, the lower bound of the excluded region does not depend on the total decay width of the scalar, being determined only by \(\varGamma ({\varphi \rightarrow \mu \bar{\mu }})\); in particular, there is no drop of the sensitivity at \(m_{S}\simeq 1\,\text {GeV}\) observed in Fig. 1 for \(\text {Br}(\varphi \rightarrow \mu \bar{\mu })\). The reason for this is that scalars with \(\sin \theta \) close to the lower bound are very long-lived, with the characteristic decay length \(c\tau _{\varphi }\langle \beta \gamma _{\varphi }\rangle \) exceeding much the geometric size of the detector. As a result, the decay probability in (5) behaves as \(1- P(r_\textrm{det}|\beta \gamma ) \propto \tau _{\varphi }^{-1}\), and the product \(\text {Br}(\varphi \rightarrow \mu \bar{\mu })\tau _{\varphi }^{-1}\) is just \(\varGamma (\varphi \rightarrow \mu \bar{\mu })\). Indeed, the scalar’s proper lifetime at the vicinity of 1 GeV and \(\sin \theta \sim 10^{-4}\) is \(\tau _{\varphi }\gtrsim {\mathcal {O}}(100)\,\text {ps}\), which corresponds to the “extremely displaced region” in the LHCb search \(\tau \gg 10\text { ps}\). Because of this, in the minimal scalar model, the lower bound does not depend on the description of the hadronic decays of the scalar (remind Fig. 1).

The situation may be different if the total scalar decay width has contributions from \(\theta \) (given by the decay width \(\varGamma ^{(\theta )}_{\varphi }\)) as well as another contribution, for instance, an invisible decay width \(\varGamma _{\textrm{inv}}\):

$$\begin{aligned} \varGamma _{\varphi ,\text {tot}}=\varGamma ^{(\theta )}_{\varphi } +\varGamma _\textrm{inv}. \end{aligned}$$
(7)

If the latter is sufficiently large to be comparable with \(\varGamma ^{(\theta )}_{\varphi }\), the lifetime becomes small enough such that all scalars decay inside the detector, see Fig. 2 (right).

The second type of searches, \(B\rightarrow K+\mathrm inv\), corresponds to the scenario when \(\varphi \) is not detected. This may happen if \(\varphi \) decays into invisible particles (such as neutrinos or hypothetical feebly-interacting particles (FIPs) that leave the detector) or if it is too long-lived and escapes the detector before decaying. Therefore, the number of events scales as

$$\begin{aligned}&\textrm{BR}(B^+\rightarrow K^+ \varphi ) \nonumber \\&\quad \times \left[ \textrm{BR}(\varphi \rightarrow \textrm{inv}) + P(r_\textrm{det}|\beta \gamma ) \, \textrm{BR}(\varphi \rightarrow f\bar{f})\right] . \end{aligned}$$
(8)

BaBar, Belle, and Belle II already searched for \(B\rightarrow K+\textrm{inv}\) [24,25,26, 42]. In the minimal scalar model all scalar decays into, e.g., \(e\bar{e},\mu \bar{\mu },\pi \pi , KK\), etc., would induce a visible activity in the detector. Therefore, the first summand is effectively zero. As for the second contribution, due to the scaling \(\tau _{\varphi }\propto \theta ^{-2}\), to be long-lived enough to escape the detector, scalars need to have small \(\theta \); otherwise, the probability is exponentially suppressed. The latter means the suppression of the production branching ratio, making “invisible” events with scalars very rare. Therefore, the second type of search is not very efficient.

However, the situation may drastically change if the scalar may decay into invisible particles, such that the invisible decay width \(\varGamma _\textrm{inv}\) becomes comparable with the widths into visible decay states. Similarly to the decays into visible particles, the missing invariant mass distribution would be peaked at \(m_{\varphi }\). This is especially useful since, in addition to the constraint on the branching ratio \(\text {Br}(B\rightarrow K+\mathrm inv)\), BaBar [24] provides the distribution in the missing squared invariant mass \(q^2\). For 2-body decays that feature a narrow resonance in the \(q^2\) distribution, almost all events are contained within one bin. This is illustrated in Fig. 3, which shows the number of events in each bin from a scalar with mixing angle \(\sin \theta =6\times 10^{-3}\) and invisible decay width \(\varGamma _\textrm{inv}=10\) eV for different scalar masses \(m_\varphi \). This will be used to extract constraints on the scalar mixing angle \(\theta \) in Sect. 4.

Fig. 3
figure 3

The plot shows the experimental BaBar data for \(B^+\rightarrow K^+ +\textrm{inv}\) [24] binned in the missing invariant mass squared normalised to the decaying B meson mass, \(s_B=q^2/m_B^2\). The background is shown as black dashed line and the SM prediction for \(B^+ \rightarrow K^+\nu \bar{\nu }\) as a solid blue line. For illustration, we indicate how a light scalar with \(\sin \theta =6\times 10^{-3}\) and \(\varGamma _\textrm{inv}=10\) eV would contribute for different masses \(m_\varphi \) using the partial decay width to SM particles calculated following [46]

4 Complementarity of \(B\rightarrow K \mu \bar{\mu }\) and \(B \rightarrow K +\textrm{inv}\)

As argued in the previous section, the simple picture changes for scalars with a sizeable invisible decay width. Decays to invisible final states leave a signal in \(B\rightarrow K +\textrm{inv}\). We recast the BaBar search for \(B^+ \rightarrow K^+ +\textrm{inv}\) [24], which is the most sensitive search channel of the different \(B\rightarrow K +\textrm{inv}\) decays. It provides the data in ten equally spaced bins in the \(q^2\) distribution as reproduced in Fig. 3 from Fig. 5 in [24]. The BaBar data is statistics limited with systematic errors at the percent level. We validated our analysis by reproducing the rescaled SM prediction and the upper bound for BR(\(B^+\rightarrow K^+ \nu \bar{\nu }\)) of the BaBar analysis [24] within 20% (19% including systematic errors).Footnote 4 As the systematic errors are negligible compared to the statistical errors and the precision of the final result, we neglect them in the following statistical analysis. Following [24] we assume the events in each bin are distributed following a Poisson distribution \(\textrm{Poisson} (k|\lambda )=\lambda ^k e^{-\lambda }/k!\) and thus the likelihood is

$$\begin{aligned} {\mathcal {L}} = \prod _{i} \textrm{Poisson}\left( N_i\Big | \epsilon _i s_i N_{B\bar{B}}+ b_i \right) \end{aligned}$$
(9)

with the signal efficiency \(\epsilon _i~\sim {\mathcal {O}}(10^{-3})\), extracted from Fig. 6 in [24], the total number of \(B\bar{B}\) events \(N_{B\bar{B}}=4.71\times 10^8\), and the expected background events in each bin, \(b_i\). The estimates for the background events are separated into peaking background events with a correctly reconstructed tagged event, estimated from Monte Carlo simulations, and combinatorial background from continuum events and incorrectly reconstructed events, which has been extrapolated from data [24]. Finally, the branching ratio for signal events in each bin is given by

$$\begin{aligned} s_i&= \tau _B \int _\mathrm{bin \,i} \!\!\!\!dq^2 \left( \frac{d\varGamma _\textrm{SM}(B^+\rightarrow K^+\nu \bar{\nu })}{dq^2} \right. \nonumber \\&\quad \left. +\frac{d\varGamma (B^+\rightarrow K^+ \varphi (\rightarrow \textrm{inv}))}{dq^2}\right) . \end{aligned}$$
(10)

The first term denotes the SM contribution, see [27, 32, 35]. As the decay width of the scalar is much smaller than its mass for the relevant parameter space, we employ the narrow width approximationFootnote 5 and find

$$\begin{aligned}&\frac{d\varGamma (B^+\rightarrow K^+ \varphi (\rightarrow \textrm{inv}))}{dq^2} =\delta (q^2-m_\varphi ^2)\varGamma (B\rightarrow K^{+} \varphi )\nonumber \\&\qquad \times \left( \frac{\varGamma _\textrm{inv}}{\varGamma _{\varphi ,\mathrm tot}} + P(r_\textrm{det}|\beta \gamma ) \frac{\varGamma _\varphi ^{(\theta )}}{\varGamma _{\varphi ,\mathrm tot}}\right) . \end{aligned}$$
(11)

The detector size is set to \(r_\textrm{det}=0.5\) m following [51]. For large invisible decay width \(\varGamma _\textrm{inv}\gg \varGamma _{\varphi }^{(\theta )}\), the second term is negligible and the differential decay width is proportional to \(\sin ^2\theta \), while for small invisible decay width the total decay width \(\varGamma _{\varphi ,\mathrm tot} \approx \varGamma _{\varphi }^{(\theta )} \propto \sin ^2\theta \), and thus the \(\sin ^2\theta \) dependence cancels in the differential decay width. Hence, the constraint from \(B^+\rightarrow K^++\textrm{inv}\) can be interpreted as a constraint on \(\varGamma _\textrm{inv}\) for small invisible decay widths and on \(\sin \theta \) for large invisible decay widths. Without performing any statistical analysis, from the few number of observed events, the total number of \(B\bar{B}\) mesons and the efficiency, we expect to be sensitive to branching ratios BR\((B^+\rightarrow K^+\varphi )\sim {\mathcal {O}}(10^{-5})\).

From the likelihood function, the corresponding \(\chi ^2\) function is given by

$$\begin{aligned} \chi ^2 = -2\ln {\mathcal {L}} = \sum _i f(N_i|\epsilon _i s_i N_{B\bar{B}} +b_i^\textrm{peak} + b_i^\textrm{comb}) \end{aligned}$$
(12)

with \(f(n|\nu )=2\nu -2 n\ln \nu +2 \ln (n!)\). Minimising the \(\chi ^2\) function with respect to the scalar mixing angle \(\sin \theta \) for fixed scalar mass \(m_\varphi \) and invisible decay width \(\varGamma _\textrm{inv}\), we derive an upper limit on \(\sin \theta \) at 90% CL (\(\varDelta \chi ^2=2.7\)).

Fig. 4
figure 4

Exclusion contours for different scalar masses as a function of the invisible decay width \(\varGamma _\textrm{inv}\). Mixing angles above the dashed lines are excluded by \(B^+\rightarrow K^+ \varphi (\rightarrow \mu \bar{\mu })\) from LHCb and above the solid lines by \(B^+\rightarrow K^+ +\textrm{inv}\) from BaBar. We use the prediction for the scalar branching ratios from [46] with included NLO corrections. The blue dotted line shows the BaBar constraint based on [48]. The double-dot-dashed (dot-dashed) orange line for \(m_\varphi =2.2\) GeV corresponds to BR(\(B^+ \rightarrow K^+ \varphi \)) being equal to (10% of) the SM prediction. The grey line indicates \(\varGamma _\textrm{inv}=\varGamma _{\varphi }^{(\theta )}\) for \(m_\varphi =2.2\) GeV. The stars indicate the maximum invisible decay width \(\varGamma _\textrm{inv}=\varGamma (\varphi \rightarrow NN)\) in the gauged \(B-L\) model

Fig. 5
figure 5

Recast LHCb exclusions (shaded regions) for a Higgs-like scalar for selected fixed invisible decay widths in comparison to the upper bounds on the mixing angle \(\theta \) from the BaBar search for \(B^+\rightarrow K^+ +\textrm{inv}\) for the same invisible decay widths using solid lines with the same colours. The partial decay width to SM particles in the left (right) plot has been calculated following [46] ([48])

The main results are presented in Figs. 4 and 5, which show the constraints on the scalar mixing angle \(\sin \theta \) from the search for \(B^+\rightarrow K^+ \varphi (\rightarrow \mu \bar{\mu })\) at LHCb [20] and invisible decay searches at BaBar [24] as a function of the invisible decay width \(\varGamma _\textrm{inv}\) and the scalar mass \(m_\varphi \) for different benchmark values. The BaBar constraints are shown as solid contours in both figures, while the LHCb constraints are indicated by dashed lines in Fig. 4 and by shaded regions in Fig. 5. For small invisible decay widths \(B^+\rightarrow K^+ \varphi (\rightarrow \mu \bar{\mu })\) places the most stringent constraint on \(\sin \theta \), while \(B^+\rightarrow K^+ \varphi (\rightarrow \textrm{inv})\) places the most stringent constraint for large invisible decay widths. The cross over occurs for \(\varGamma _\textrm{inv} \sim {\mathcal {O}}(1{-}50)\) eV depending on the scalar mass.

The dependence of the LHCb search for \(B^+\rightarrow K^+\varphi (\rightarrow \mu \bar{\mu })\) [20] is straightforward to understand. For \(\varGamma _\textrm{inv}\ll \varGamma _{\varphi }^{(\theta )}\), there is no dependence on the invisible decay width, as can be observed on the left side of Fig. 4. The grey line shows \(\varGamma _\textrm{inv} =\varGamma _\varphi ^{(\theta )}\) for \(m_\varphi =2.2\) GeV (orange contours). At the intersection of the dashed orange and grey line, the LHCb constraint starts to weaken. For \(\varGamma _\textrm{inv} \gtrsim 3\) meV, the LHCb constraints are described by lines on a log-log scale in Fig. 4, because the scalar decays promptly, while for \(\varGamma _\textrm{inv}\lesssim 3\) meV, the constraint depends on the finite scalar lifetime. Similarly in Fig. 5, the constraint for \(\varGamma _\textrm{inv}=1\) keV can be obtained from the one for \(\varGamma _\textrm{inv}=1\) eV by rescaling, while the constraints for \(\varGamma _\textrm{inv}\lesssim 10\) meV feature a non-trivial dependence on the invisible decay width \(\varGamma _\textrm{inv}\) via the dependence on the scalar lifetime \(\tau _{\varphi }\). Note, the LHCb result is largely insensitive to the different predictions [46, 48] for the branching ratios.

The dependence of the BaBar constraints (solid contours) in Fig. 4 can be understood from the above argument: For large \(\varGamma _\textrm{inv}\gg \varGamma _{\varphi }^{(\theta )}\), \(B\rightarrow K+\textrm{inv}\) can be interpreted as a constraint on \(\sin \theta \), while for small \(\varGamma _\textrm{inv}\) it has to be interpreted as a constraint on \(\varGamma _\textrm{inv}\), which explains the sharp drop in sensitivity for \(\varGamma _\textrm{inv}\lesssim 5\) meV. This can also be observed in Fig. 5: The BaBar constraints for \(\varGamma _\textrm{inv}=1\) eV and 1 keV agree except for scalar masses close to the kinematic cutoff. There is no sensitivity for \(\varGamma _\textrm{inv}=0\) eV and we observe a strong dependence on the scalar mass for \(\varGamma _\textrm{inv}=10\) meV, because the BaBar sensitivity depends on the \(q^2\) bin. To illustrate the sensitivity of \(B^+\rightarrow K^+ +\textrm{inv}\) to new physics, we show iso-contours for the branching ratio of BR(\(B^+\rightarrow K^+\varphi )\) being equal to (10% of) the SM branching ratio BR(\(B^+\rightarrow K^+ \nu \bar{\nu }\)) for a scalar with \(m_{\varphi }=2.2\) GeV as dot-dashed (dot-dot-dashed) orange lines in Fig. 4.

Finally, depending on the scalar mass, the theoretical uncertainty in the scalar’s hadronic decay width affects the lower bounds for visible/invisible signatures. Let us first consider the invisible case. In Fig. 4, we consider two descriptions of the scalar’s width discussed in Sect. 2 for the mass \(m_{\varphi } = 0.9\) GeV: as in [46] (the solid blue line) and [48] (the dotted blue line). Also, Fig. 5 shows the sensitivity from the invisible signature assuming these two descriptions (the left and right panels correspondingly). In the domain of large \(\varGamma _\textrm{inv}\gg \varGamma ^{(\theta )}_{\varphi }\), where \(\varGamma ^{(\theta )}_{\varphi }\) is the total width controlled by \(\theta \), the dependence of the number of events on the total width disappears. Therefore, the sensitivity does not depend on the uncertainty, as we see for \(\varGamma _\textrm{inv}\). Indeed, in Eq. (8), the first summand in the brackets reduces to 1, while the second summand tends to 0. Once \(\varGamma _\textrm{inv}\) decreases, \(\varGamma ^{(\theta )}_{\varphi }\) becomes essential. In particular, in the limit \(\varGamma _\textrm{inv}\ll \varGamma ^{(\theta )}_{\varphi }\), the number of events scales as \((\varGamma ^{(\theta )}_{\varphi })^{-1}\): larger width means lower number of events. The width from [46] is resonantly enhanced compared to the one from [48], and therefore the sensitivity of the invisible signature is weaker in the former case.

For the visible case, the number of events scales as given in Eq. (5). We illustrate the results for two different descriptions in Fig. 5. As we already discussed, if the invisible width is very small or zero, the LHCb bounds extend to small values of the mixing angle where the probability of the scalar to decay inside the decay cancels with the \(\text {Br}(\varphi \rightarrow \mu \bar{\mu })\). In the opposite case of large \(\varGamma _\textrm{inv}\), the decay probability does not have such scaling and tends to 1. The dependence on the total width via \(\text {Br}(\varphi \rightarrow \mu \bar{\mu })\) survives. Similarly to the invisible signature, larger \(\varGamma ^{(\theta )}_{\varphi }\) means weaker bounds. In particular, assuming the description from [46] and \(\varGamma _\textrm{inv} = 1\,\text {eV}\), we see that the sensitivity is reduced at \(m_\varphi =1\) GeV due to the resonant enhancement of the \(f_{0}(980)\) and thus large visible decay width into SM particles via the scalar mixing. The shape is different if assuming the description as in [48], where no resonant enhancement exists.

5 SM extensions with a light scalar with invisible width

There are several possibilities how the light scalar may decay and escape detection in the LHCb searches for \(B\rightarrow K \varphi (\rightarrow \mu \bar{\mu })\). An attractive possibility is to couple the scalar field to fermionic dark matter. However thermal production of this dark matter candidate in the early universe is strongly constrained [2] and production via freeze-in requires small couplings. We thus focus on scalar decays to unstable particles which further decay to lighter SM particles, like HNLs. HNLs are well-motivated as explanation of neutrino masses via the seesaw mechanism [52]. Big bang nucleosynthesis (BBN) places a lower bound on the lifetime of GeV-scale HNLs of \(\tau _N<0.02\)s [53] and thus a lower bound on the active-sterile mixing angle. Direct searches on the other hand provide an upper bound on the active-sterile mixing angle. Assuming that the HNLs generate neutrino masses, together the two constraints exclude HNL masses below 0.33–0.36 GeV, set by the kinematic threshold of \(K\rightarrow \pi \nu N\), apart from a small window \(M_N\in [0.12,0.14]\) GeV [54]. For the allowed parameter space GeV-scale HNLs escape the detector undetected, see e.g. [54, 55].

5.1 Real scalar coupling to heavy neutral leptons

The simplest viable scenario is to couple a real scalar field to HNLs \(N_i\), which allows to generate the invisible scalar decay width via \(\varphi \rightarrow NN\). The relevant Lagrangian is

$$\begin{aligned} {\mathcal {L}} =- \frac{1}{2} N^TC (\mu _N + y_N \varphi ) N + \mathrm{h.c.} \end{aligned}$$
(13)

where C is the charge conjugation matrix, \(M_N = \mu _N+ y_N v_\phi \) is the HNL mass term and \(y_N\) the Yukawa coupling to the scalar \(\varphi \). Hence, the Yukawa coupling \(y_N\) and the HNL mass can be chosen independently.

An important constraint on this scenario comes from SM Higgs decays. In presence of a quartic Higgs portal interaction \(\tfrac{\lambda _{H\varphi }}{2} H^\dagger H \varphi ^2\), the SM Higgs can decay in the real scalar with the branching ratio

$$\begin{aligned} \textrm{BR}(h\rightarrow \varphi \varphi ) = \frac{\lambda ^{1/2} (m_h^2,m_\varphi ^2,m_\varphi ^2)}{32 \pi \sqrt{2} G_F m_h^3 \varGamma _h} |\lambda _{H\varphi }|^2 \end{aligned}$$
(14)

with the Källén function \(\lambda (x,y,z) =x^2+y^2+z^2-2xy-2xz-2yz\). The upper bound BR(\(h\rightarrow \textrm{inv})\le 0.18\) [56] constrains the quartic coupling as \(|\lambda _{H\varphi }|\lesssim 0.01\). Similarly, for non-zero vacuum expectation value (VEV) of \(\varphi \), the quartic Higgs portal interaction contributes to the scalar mixing. Demanding this contribution to be small results in \(\lambda _{H\varphi }\ll \sqrt{2} \sin \theta m_h^2/vv_\varphi \) with the electroweak VEV \(v^2 =1/\sqrt{2} G_F\), i.e. we assume that the mixing is dominated by the cubic term in the Higgs potential.

Even in absence of a quartic Higgs portal interaction, there is a contribution to invisible Higgs decay from \(h\rightarrow NN\) which translates into an upper bound on the HNL Yukawa coupling

$$\begin{aligned} \textrm{BR}(h\rightarrow N N)&= \sum _{i,j}\frac{|y_{Nij}|^2 \sin ^2 \theta }{8\pi (1+\delta _{ij})} \frac{m_h}{\varGamma _h} \nonumber \\&\Rightarrow \sum _{i,j} \frac{|y_{Nij}|^2}{1+\delta _{ij}} \lesssim 1.4 \,\left( \frac{10^{-2}}{\sin \theta }\right) ^2. \end{aligned}$$
(15)

The partial decay width for scalar decay to HNLs is

$$\begin{aligned} \varGamma (\varphi \rightarrow N_iN_j)&= \frac{m_\varphi | y_{Nij}|^2 }{8\pi (1+\delta _{ij})} \lambda ^{1/2}(1,x_i^2,x_j^2)\nonumber \\&\quad \times (1-(x_i+x_j)^2) \end{aligned}$$
(16)

with \(x_i\equiv M_{Ni}/m_\varphi \). Taking the constraint from invisible Higgs decay into account, we obtain an upper bound for the invisible width for \(x_i\ll 1\)

$$\begin{aligned} \varGamma (\varphi \rightarrow NN)&\lesssim 0.06 \,\left( \frac{10^{-2}}{\sin \theta }\right) ^2\, m_\varphi . \end{aligned}$$
(17)

Hence, there is no strong constraint on the partial decay width \(\varGamma (\varphi \rightarrow NN)\) and a sizeable invisible decay width is allowed for a real scalar coupling to HNLs. As argued above, the lifetime of HNLs is long compared to the size of the detector and thus they escape undetected. For example, using the lifetime of HNLs from [57], we explicitly find for \(m_N=2.3\) GeV a decay length of \(c\tau _N=3.75\times 10^5\) m at the seesaw line \(U^{2}_{\textrm{seesaw}} \sim 5\times 10^{-11} (1\text { GeV}/m_{N})\) [58]. Together with the allowed mass range of HNLs, this demonstrates the whole range of invisible partial decay widths in Fig. 4 can be obtained in the real singlet model. The model is also able to explain light neutrino masses via the seesaw mechanism [52].

5.2 \(B-L\) gauge symmetry

The real scalar field discussed in the previous subsection has many undetermined parameters, which can be reduced by introducing a symmetry. A well-motivated scenario is the \(B-L\) symmetry. After spontaneous breaking of the \(B-L\) symmetry, a Majorana mass term for the HNLs is generated, which provides an explanation of active neutrino masses via the seesaw mechanism [52]. A global \(B-L\) symmetry is straightforwardly ruled out: spontaneous breaking of the \(B-L\) symmetry results in a Majoron, which is efficiently produced via its interactions with the HNLs. It thus contributes to the effective number of neutrinos, \(N_\textrm{eff}\), and as it only decouples below the scale of the HNLs, its contribution to \(N_\textrm{eff}\) is too large and excluded by BBN. These constraints are avoided by promoting the global \(B-L\) symmetry to a gauge symmetry, which has been first proposed in [10, 11].

The relevant interactions for the following discussion are

$$\begin{aligned} {\mathcal {L}}&\supset (D_\mu \phi )^\dagger ( D^\mu \phi ) - \lambda _{H\phi }\left( |\phi |^2-\frac{v_\phi ^2}{2}\right) \left( H^\dagger H -\frac{v^2}{2}\right) \nonumber \\&\quad -\frac{1}{2} N^TC y_N \phi N. \end{aligned}$$
(18)

After spontaneous breaking of the \(B-L\) gauge symmetry with \(\langle \phi \rangle = v_\phi /\sqrt{2}\), HNLs and the \(Z'\) gauge boson acquire the massesFootnote 6\(M_{Ni}= y_{Ni} v_\phi /\sqrt{2}\) and \(m_{Z'}= 2 g_\textrm{BL} v_\phi \) with the \(B-L\) gauge coupling. The interactions of the real scalar \(\varphi \equiv \textrm{Re}(\phi )/\sqrt{2}\) are by construction proportional to the HNL and \(Z'\) masses

$$\begin{aligned} {\mathcal {L}} \supset - \frac{1}{2} \varphi N_i^T C\frac{M_{Ni}}{v_\phi } N_i + \frac{1}{2} \frac{m_{Z'}^2}{v_\phi ^2}\varphi ^2 Z^\prime _\mu Z^{\prime \mu } + \frac{m_{Z'}^2}{v_\phi }\varphi Z^\prime _\mu Z^{\prime \mu } \end{aligned}$$
(19)

which results in tight connections between the different observables.

The partial decay width for decay to HNLs in Eq. (16) can be expressed in terms of the masses and the gauge coupling \(g_\textrm{BL}\)

$$\begin{aligned} \varGamma (\varphi \rightarrow NN) = \frac{g_\textrm{BL}^2 m_\varphi ^3}{2\pi m_{Z'}^2} \sum _i x_i^2 (1-4x_i^2)^{3/2} \end{aligned}$$
(20)

and the decay to a pair of \(Z'\) gauge bosons is given in terms of

$$\begin{aligned} \varGamma (\varphi \rightarrow Z'Z') = \frac{g_\textrm{BL}^2 m_\varphi }{8\pi } \frac{(1-4z)^{1/2}\left( 1 -4z +12z^2 \right) }{z} \end{aligned}$$
(21)

with \(z=m_{Z'}^2/m_\varphi ^2\). The \(Z'\) gauge boson however is not stable and further decays to SM particles. Its lifetime is inversely proportional to the square of the gauge coupling, \(\tau _{Z'} \propto g_\textrm{BL}^{-2} m_{Z'}^{-1}\). Hence, the \(Z'\) gauge boson quickly decays to neutrinos, and depending on its mass to charged leptons and hadrons. The standard model charged leptons couple to \(Z'\) vectorially and the branching ratio for active neutrinos can be straightforwardly obtained as

$$\begin{aligned} \textrm{Br}(Z' \rightarrow \ell \bar{\ell })&=\frac{g_\textrm{BL}^2}{12\pi } \frac{m_{Z'}}{\varGamma _{Z'}}\left( 1-4y\right) ^{1/2} \left( 1+2 y\right) \nonumber \\ \textrm{Br}\left( Z' \rightarrow \sum _i\nu _i \bar{\nu }_i \right)&=\frac{g_\textrm{BL}^2}{8\pi }\frac{m_{Z'}}{\varGamma _{Z'}} \end{aligned}$$
(22)

with \(y=m_\ell ^2/m_{Z'}^2\). Hence, a substantial fraction of the \(Z'\) gauge bosons decays to visible final states which leave a signal in the detector with multiple leptons and/or hadrons and thus do not contribute to either \(B^+\rightarrow K^+\varphi (\rightarrow \mu \bar{\mu })\) or \(B^+\rightarrow K^+ \varphi (\rightarrow \textrm{inv})\). Variants of the \(B-L\) gauge symmetry with non-universal lepton number, see e.g. [59, 60], may forbid \(Z'\) decays to light leptons and thus evade any constraints from the cascade decay \(B^+\rightarrow K^+ \varphi (\rightarrow Z' Z'\rightarrow \ell \bar{\ell } \ell \bar{\ell })\), see e.g. [61] for a discussion of a decay with two muon–antimuon pairs in the final state. Final states with multiple muons have also been studied in [62]. In the following we consider the scenario, where the decay \(\varphi \rightarrow Z'Z'\) is kinematically forbidden by choosing \(m_{Z'}\ge m_\varphi /2\).

Fig. 6
figure 6

Maximum decay width \(\varGamma (\varphi \rightarrow NN)\) as function of \(m_\varphi \) and \(m_{Z'}/g_\textrm{BL}\). The grey-shaded region is excluded by \(Z'\) searches assuming \(m_{Z'}\ge m_\varphi /2\) such that \(\varphi \rightarrow Z'Z'\) is kinematically forbidden

In Fig. 6 we present the maximum decay width \(\varGamma (\varphi \rightarrow NN)\) as a function of the scalar mass \(m_\varphi \) and the \(Z'\) mass \(m_{Z'}/g_\textrm{BL}\). For this we consider two degenerate HNLs with masses \(m_N=\max (\tfrac{m_\phi }{\sqrt{10}},0.3\, \textrm{GeV})\) which maximises the partial decay width \(\varGamma (\varphi \rightarrow NN)\) and respects the lower bound on the mass of HNLs. The grey-shaded region is excluded by searches for light \(Z'\) bosons. It has been extracted from the top-left plot of Fig. 2 in [63] (see also [64, 65]). The coloured solid contours show \(\varGamma _\textrm{inv} \in [0.01, 0.1,1,10]\) eV. Together this demonstrates the possibility of a sizeable invisible scalar decay width \(\varGamma _\textrm{inv}\gtrsim 1\) eV. We indicate the maximum possible invisible decay width consistent with the constraint on \(m_\mathrm{Z'}/g_\textrm{BL}\) also in Fig. 4 as stars for the three benchmark masses.

Finally, invisible Higgs decays place a constraint on the gauge \(B-L\) model via the quartic Higgs portal \(H^\dagger H \phi ^\dagger \phi \). For small mixing and light scalars, \(m_\varphi \ll m_h\), the Higgs branching ratio to a pair of scalars is given by

$$\begin{aligned} \textrm{BR}(h\rightarrow \varphi \varphi ) = \frac{ g_\textrm{BL}^2}{8\pi } \frac{ m_h^3 }{m_{Z'}^2 \varGamma _h} \sin ^2\theta , \end{aligned}$$
(23)

where we expressed the Higgs portal coupling in terms of the scalar mixing angle \(\theta \). Higgs decays to HNLs and \(Z'\) gauge bosons are negligible in this case, because they are suppressed by the small masses of the final state particles. The constraint on the invisible Higgs decay branching ratio thus results in an upper bound on \(\sin \theta \),

$$\begin{aligned} \sin ^2\theta \le 0.93 \left( \frac{m_{Z'}}{\textrm{GeV}} \right) ^2 \left( \frac{10^{-4}}{g_\textrm{BL}}\right) ^2. \end{aligned}$$
(24)

It is however weaker than the searches for the scalar in B meson decays and does not pose any additional constraint.

6 Conclusions

We considered a scalar, which mixes with the SM Higgs boson with an additional invisible decay width \(\varGamma _\textrm{inv}\), and scrutinised the constraints on the scalar mixing angle from the search \(B^+\rightarrow K^+ \varphi (\rightarrow \mu \bar{\mu })\) at LHCb [20]. The LHCb search looses sensitivity for invisible decay widths larger than the visible decay width, \(\varGamma _\textrm{inv} \gtrsim \varGamma _{\varphi }^{(\theta )}\). We demonstrated that for scalars with an invisible decay width of \(\varGamma _\textrm{inv}\gtrsim {\mathcal {O}} (1{-}50)\) eV, the decay \(B^+\rightarrow K^+ \varphi (\rightarrow \textrm{inv})\) provides the most stringent constraint on the scalar mixing angle \(\theta \), which opens an opportunity for Belle II to discover new physics in \(B^+\rightarrow K^+ +\textrm{inv}\).

We provided two explicit models which realise a sizeable invisible decay width. The gauged \(B-L\) model is well motivated as an explanation of non-zero neutrino masses. In this model, the scalar breaking \(B-L\) gauge symmetry may decay to heavy neutral leptons. The heavy neutral leptons are long-lived on the time scale of the detector, escape it undetected and thus contribute to the invisible decay width of the \(B-L\) scalar. This scenario is mostly constrained by searches for the \(Z'\) gauge boson which limits the invisible decay width to less than \({\mathcal {O}}(30)\) eV. As those constraints do not apply in the real scalar model, it is possible to obtain an invisible decay width in the MeV range.

Finally, as the phenomenology mainly depends on scalar mixing and the coupling of the scalar to heavy neutral leptons, the conclusions from the \(B-L\) model apply at least qualitatively to other gauged \(U(1)'\) extensions of the SM, as long as the scalar spontaneously breaking the \(U(1)'\) symmetry can decay to heavy neutral leptons.