BBN for the LHC: constraints on lifetimes of the Higgs portal scalars

TL;DR

This paper investigates the lifetime constraints of Higgs portal scalars from Big Bang Nucleosynthesis, showing most must decay within 0.1 seconds to avoid cosmological conflicts, impacting LHC and MATHUSLA searches.

Contribution

It provides the first detailed analysis of BBN constraints on the lifetimes of Higgs portal scalars within a minimal model, linking cosmology with collider phenomenology.

Findings

Lifetimes must be less than 0.1 seconds for scalars with mass between 2m_mu and half the Higgs mass.

Lifetimes can extend to about one second for scalars below 2m_mu.

Strong BBN constraints support the feasibility of detecting displaced decays at MATHUSLA.

Abstract

LHC experiments can provide a remarkable sensitivity to exotic metastable massive particles, decaying with significant displacement from the interaction point. The best sensitivity is achieved to models where the production and decay occur due to different coupling constants, and the lifetime of exotic particles determines the probability of decay within a detector. The lifetimes of such particles can be independently limited from standard cosmology, in particular the Big Bang Nucleosynthesis. In this paper, we analyze the constraints on the simplest scalar model coupled through the Higgs portal, where the production occurs via $h\to SS$, and the decay is induced by the small mixing angle of the Higgs field $h$ and scalar $S$. We find that throughout the most part of the parameter space, $2 m_\mu < m_S < m_h/2$, the lifetimes of exotic particle has to be less than 0.1 seconds, while below $2m_\mu$ it could grow to about a second. The strong constraints on lifetimes are induced by the nucleonic and mesonic decays of scalars that tend to raise the $n/p$ ratio. Strong constraints on lifetimes of the minimal singlet extensions of the Higgs potential is a welcome news for the MATHUSLA proposal that seeks to detect displaced decays of exotic particles produced in the LHC collisions. We also point out how more complicated exotic sectors could evade the BBN lifetime constraints.

Equations

$${\mathcal{L}}_{H/S}={\mu }^2{H}^{†}H-{\lambda }_H{\left(H^{†}H\right)}^2-V(S)-AS{H}^{†}H-{\lambda }_S{S}^2{H}^{†}H+\text{mbox}kin.terms.$$

$$\theta =\frac{Av}{m_h^2-m_S^2}\left(1-\frac{{\lambda }_S{v}^2}{m_S^2}\right).$$

$${\mathcal{L}}_{\mathrm{decay}}=S\times \theta \sum _{\mathrm{SM}}{O}_h,$$

$${\Gamma }_{h\to SS}$$

$$Br(h\to SS)$$

$${\lambda }_S\lesssim \frac{0.007}{{\left(1-\frac{4m_S^2}{m_h^2}\right)}^{1/4}}.$$

$${\Gamma }_{S\to l\bar{l}}=\frac{{\theta }^2}{8\pi }\frac{m_l^2}{v^2}{m}_S{\left(1-\frac{4m_l^2}{m_S^2}\right)}^{3/2}.$$

$$f_q=3\left[1+\frac{4{\alpha }_s}{3\pi }\left(\frac{9}{4}+\frac{3}{2}\ln \frac{m_q^2}{m_S^2}\right)\right]$$

$${\mathcal{L}}_{S\pi \pi }=\frac{4}{9}\frac{\theta }{v}S\left(\frac{1}{2}{\partial }_{\mu }{\pi }^0{\partial }^{\mu }{\pi }^0+{\partial }_{\mu }{\pi }^{+}{\partial }^{\mu }{\pi }^{-}\right)-\frac{5}{3}\frac{\theta m_{\pi }^2}{v}S\left(\frac{1}{2}{\pi }^0{\pi }^0+{\pi }^{+}{\pi }^{-}\right),$$

$${\Gamma }_{S\to {\pi }^{+}{\pi }^{-}}=2{\Gamma }_{S\to {\pi }^0{\pi }^0}=\frac{{\theta }^2}{16\pi }\frac{m_S^3}{v^2}{\left(\frac{2}{9}+\frac{11}{9}\frac{m_{\pi }^2}{m_S^2}\right)}^2\sqrt{1-\frac{4m_{\pi }^2}{m_S^2}}.$$

$${\Gamma }_{{\mu }^{+}{\mu }^{-}}:{\Gamma }_{KK}:{\Gamma }_{\eta \eta }=m_{\mu }^2{\beta }_{\mu }^3:3\frac{9}{13}{m}_s^2{\beta }_K^3:3\frac{4}{13}{m}_s^2{\beta }_{\eta }^3,$$

$$\sigma v(s)=\frac{8{\lambda }_S^2{v}^2}{(s-m_h^2{)}^2+m_h^2{\Gamma }_{\mathrm{SM}+S}^2}\frac{{\Gamma }_{\mathrm{SM}}^{m_h\to \sqrt{s}}}{\sqrt{s}},⟨\sigma v⟩=\frac{{\int }_{4m_S^2}^{\infty }ds\sigma v(s)s\sqrt{s-4m_S^2}{K}_1\left(\frac{\sqrt{s}}{T}\right)}{16T{m}_S^4{K}_2^2\left(\frac{m_S}{T}\right)}.$$

$$\frac{dY}{dx}=\frac{s⟨\sigma v⟩}{Hx}\left[1+\frac{1}{3}\frac{d(\ln h_{\mathrm{eff}})}{d(\ln T)}\right]\left(Y_{\mathrm{eq}}^2-Y^2\right),$$

$$P_{n\to p}\simeq n_p⟨\sigma v{⟩}_{pn}c{\tau }_{{\pi }^{+}}\simeq 2\times \frac{10^{21}}{\mathrm{c}{\mathrm{m}}^3}\times 1.5\mathrm{mb}\times 2.6\times 10^{-8}\sec \times c\simeq 2.5\times 10^{-3}.$$

$${\Gamma }_{n{\nu }_e\to p{e}^{-}}$$

$${\Gamma }_{n{e}^{+}\to p{\bar{\nu }}_e}$$

$$J(a,b)$$

$$\frac{d{X}_n}{dT}=\frac{{\Gamma }_{n{\nu }_e\to p{e}^{-}}+{\Gamma }_{n{e}^{+}\to p{\bar{\nu }}_e}}{TH(T)}\left(X_n-(1-X_n)e^{-Q/T}\right)+\frac{{\Gamma }_n{X}_n}{TH(T)},$$

$$\Delta Y_p\equiv |Y_p-Y_p^{\mathrm{SBBN}}|<0.01,$$

$${\pi }^{-}+p$$

$${\pi }^{+}+n$$

$${\frac{d{X}_n}{dT}|}_{\pi -}=\frac{-1}{TH(T)}{n}_{{\pi }^{-}}^{\mathrm{inj}}\left(⟨\sigma v{⟩}_{pn({\pi }^0)}^{{\pi }^{-}}+⟨\sigma v{⟩}_{pn(\gamma )}^{{\pi }^{-}}\right)(1-X_n),$$

$${\bar{K}}^0+n\to K^{-}+p:(\sigma v{)}_{pn(K^{-})}^{{\bar{K}}^0}\simeq 10\mathrm{mb},Q=5.3\mathrm{MeV},$$

$$K^{-}+p$$

$$K^{-}+n$$

$$K_L+p$$

$$K_L+n$$

$$P_{p\bar{p}}^{n\to p}=\frac{X_n⟨\sigma v{⟩}_{n\bar{p}}}{X_n⟨\sigma v{⟩}_{n\bar{p}}+(1-X_n)⟨\sigma v{⟩}_{p\bar{p}}},P_{p\bar{n}}^{n\to p}=\frac{(1-X_n)⟨\sigma v{⟩}_{p\bar{n}}}{X_n⟨\sigma v{⟩}_{n\bar{n}}+(1-X_n)⟨\sigma v{⟩}_{p\bar{n}}}$$

$$⟨\sigma v{⟩}_{n\bar{n}}=⟨\sigma v{⟩}_{p\bar{p}}/C=37\mathrm{mb},⟨\sigma v{⟩}_{n\bar{p}}=⟨\sigma v{⟩}_{p\bar{n}}=28\mathrm{mb},$$

$${\frac{d{X}_n}{dT}|}_{pn}=\frac{-{\xi }_p{\Gamma }_S{e}^{-t{\Gamma }_S}}{TH(T)}\left(-P_{p\bar{p}}^{n\to p}-\kappa P_{p\bar{n}}^{n\to p}+\kappa P_{p\bar{n}}^{n\to p}+{\kappa }^2{P}_{n\bar{n}}^{n\to p}\right).$$

$${\Gamma }_{\mathrm{coll}}^{{\nu }_e}\left(E_{\nu },T\right)$$

$$\simeq \left(\frac{E_{\nu }}{32\mathrm{MeV}}\right)\left[\frac{5.7}{\sec }{\left(\frac{T_{\nu }}{1\mathrm{MeV}}\right)}^4+\frac{1.3}{\sec }\eta (T_{\gamma }){\left(\frac{T_{\gamma }}{1\mathrm{MeV}}\right)}^4\right],$$

$$f_e\left(T,E_{\nu }\right)={\Gamma }_S{Y}_S{\int }_T^{\infty }\frac{d{T}_1{e}^{-t_1{\Gamma }_S}}{H(T_1)T_1}{F}_e\left(E_{\nu },\frac{E_{0}T}{T_1}\right)e^{-{\int }_T^{T_1}d{T}_2\frac{{\Gamma }_{\mathrm{coll}}(E_{\nu }\frac{T_2}{T},T_2)}{H(T_2)T_2}},$$

Full text

I Introduction

The absence of readily discoverable new physics (NP) at the LHC has presented the physics community with a formidable puzzle. While the arguments for NP “not too far” from the weak scale still loom large, there is a distinct desire to explore wider (and wilder) theoretical options away from a simply realized weak-scale supersymmetry, or extra space dimensions. One possible strategy to look for new physics is to abandon theoretical preconceptions, and start looking for non-standard signatures that the NP could present.

Large classes of models offer promising avenues for a non-standard signal in the production of new exotic particles (possibly of electroweak-scale mass) with subsequent decay away from the interaction point (see e.g. Graham:2012th ; Craig:2015pha ; Izaguirre:2015pga ; Batell:2016zod ). While both ATLAS and CMS have performed corresponding studies in a variety of contexts and for different ranges of displacement CMS:2014hka ; Aad:2015rba ; Aad:2015uaa , it has been recently pointed out that a dedicated and relatively inexpensive detector Chou:2016lxi could extend the physics reach into cases where the decay lengths are on the order of $O(100\,{\rm m})$ and beyond.

When both the production and decay of an exotic state $S$ occur through one and the same coupling constant, the chances of detecting such NP at the LHC experiments are not great. Indeed, a large displacement implies a very small value for the coupling, which in turn leads to very inefficient production rates. Therefore, an ideal case for the collider studies would be when the production and decay occur through different coupling constants, and $\lambda_{\rm production}\gg\lambda_{\rm decay}$. For the pair-produced exotics, such a hierarchy can be made “natural” as the $\lambda_{\rm decay}\to 0$ limit could lead to an enhanced symmetry.

If the main signal to search for is an appearance of abnormal energy deposition or exotic vertex some distance from the interaction point, it is then very important to know how small $\lambda_{\rm decay}$ is allowed to be. In more practical terms, one would like to know if there is an external to the LHC physics constraint on the lifetimes $\tau_{S}$ of such exotic particles. An obvious source for such a constraint can be early cosmology. The big bang nucleosynthesis (BBN), and its overall agreement with observations Cyburt:2015mya (apart from the unclear status of 7Li) can provide a limit on the lifetimes of such particles. In order to derive such limits, one would have to make a fairly natural assumption that the Universe was indeed as hot as $T\sim m_{S}\sim~{}$electroweak scale at some point in its history. Subsequent thermal evolution to the BBN temperatures involves self-depletion via $SS\to{\rm SM}$ due to $\lambda_{\rm production}$, in an expected WIMP-type annihilation process, and late-time decay of $S\to{\rm SM}$ where depending on lifetimes and the properties of the decay products the BBN outcome may get affected. These mechanisms are well-understood in the BBN literature (see e.g. Jedamzik:2009uy ; Pospelov:2010hj for reviews). We will require that the late decay of $S$ provides a small and acceptable perturbation to the standard BBN (SBBN) outcome, which in turn will limit $\tau_{S}$.

In this paper, we analyze a fairly minimal model, where a new singlet scalar has predominantly a quadratic coupling to the Higgs boson that regulates both its production at colliders and the intermediate cosmological abundance at $T_{\rm BBN}\ll T\ll m_{S}$. Given that the model is very predictive, it allows to place robust bounds on lifetimes of such particles with a minimum amount of model dependence. We find that for most of the analyzed parameter space with $m_{S}<m_{h}/2$, the intermediate abundance of such particles is large enough to affect the neutron-proton freeze out ratios at relevant temperatures. This allows us to set fairly robust bounds on lifetimes of such particles, which come out to be remarkably strong, and shorter than $0.1$ seconds. In what follows we describe the model and the cosmological history of $S$ (section II); derive the impact on the BBN (section III); present our results (section IV), and provide related discussion (section V).

II The minimal Higgs portal model

We consider the simplest extension of the SM by a singlet scalar field $S$. A new singlet scalar $S$ can have two interaction terms with the Standard Model (SM) at the renormalizable level, in addition to trilinear and quartic self-interactions. In this scenario, the Lagrangian of the singlet sector (including the SM) generically takes the form

[TABLE]

The Higgs expectation value $v=246~{}{\rm GeV}$ is assumed to correspond to a global minimum. The self-interaction potential $V(S)=\lambda_{4}S^{4}+\lambda_{3}S^{3}+\frac{m_{S0}^{2}}{2}S^{2}$ can be redefined in such a way that the linear term is absent. It is important that the $A,~{}\lambda_{3}\to 0$ and $\langle S\rangle=0$ limit would correspond to the case of stable $S$ particles. To simplify the discussion without sacrificing much generality, we take $\lambda_{3,4}\to 0$ and assume $Av\ll m_{S0}^{2},~{}\lambda_{S}v^{2}$.

The physical mass of $S$ receives a contribution from the electroweak symmetry breaking, $m_{S}=\sqrt{m_{S0}^{2}+\lambda_{S}v^{2}}$. At linear order in $A$, the mixing angle $\theta$ between physical excitations $S$ and $h$ is

[TABLE]

The $\lambda_{S}$ term arises because the $S$ field develops a small $A$-controlled vacuum expectation value. The mixing parameter $\theta$ leads, via the $A$ coupling constant, to the decay of $S$ particles, which can be readily derived from

[TABLE]

where $O_{h}$ is the set of the standard Higgs interaction terms, with the Higgs field removed: e.g. $O_{h}=(m_{f}/v)\bar{f}f$ for an elementary SM fermion $f$.

This Yukawa-type coupling to the SM has been tested in rare meson decays OConnell:2006rsp ; Pospelov:2007mp ; Batell:2009jf ; Schmidt-Hoberg:2013hba ; Clarke:2013aya and in proton fixed-target experiments Alekhin:2015byh . The model is mostly ruled out for large mixing angles $\theta\gtrsim 10^{-4}-10^{-2}$ over the $m_{S}\sim$ MeV - 5 GeV mass range. The proposed experiment SHiP could potential improve current sensitivity down to $\theta\sim 10^{-6}$ for $m_{S}\sim$ few GeV Alekhin:2015byh .

In the limit of $\theta\to 0$, $S$ is stable and could be the dark matter Silveira:1985rk ; McDonald:1993ex ; Burgess:2000yq . Various limits arise from searches in direct and indirect detection if the particle is stable (see Refs. Cline:2013gha ; Athron:2017kgt for recent reviews), but $\lambda_{S}$ is generically bounded from the constraints on invisible Higgs decay, independently of the direct detection limits. The Standard Model Higgs has a well-predicted decay rate into SM particles of $\Gamma_{SM}=4.07$ MeV. So far, the properties of 125 GeV resonance are remarkably consistent with the SM Higgs, and therefore there is little doubt that its width is close to $\Gamma_{SM}$. The invisible branching ratio of Higgs decay to $SS$ final state is

[TABLE]

where in the last line we assumed $Br(h\to SS)\ll 1$ and $m_{S}\ll m_{h}$. The experimental upper bound on the invisible branching ratio of a SM Higgs is 0.19 (at $2\sigma$) Belanger:2013xza , which translates into an upper bound on $\lambda_{S}$

[TABLE]

If $S$ is to be stable, such small couplings would lead to an excessive abundance of $S$, which invalidates the $Z_{2}$ symmetric case, and forces us to include the decay term. From now on, we will consider $\theta\neq 0$, or in other words the case of unstable $S$ particles. Since our analysis is motivated by the LHC physics, we will use $Br(h\to SS)$ as an input parameter, and substitute $\lambda_{S}$ everywhere employing (4) and (5).

II.1 Decay products

Since $S$ interacts with the SM in the same fashion as the Higgs with an additional $\theta$ mixing factor (3), its decay properties are similar to those of a light Higgs boson. For the derivations of the actual limits on the lifetime of $S$, we need to know its mesonic and nucleonic decay branching ratios.

The decay channels of a light Higgs have been considered in the early years of the Weinberg-Salam electroweak model Ellis:1975ap , with additional refinements as SM particles, hadronic resonances were being discovered and final-state interactions better understood Ellis:1979jy ; Raby:1988qf ; Truong:1989my . Hadronic decays in the mass range $2m_{\pi}<m_{S}\lesssim 4{~{}\rm GeV}$ are still poorly understood, with models varying by as much as a few orders of magnitude near the di-kaon threshold Clarke:2013aya .

The leptonic decay channels are straightforward, with the decay rate given by

[TABLE]

If the decaying product is a pair of heavy quarks, there are $\mathcal{O}(1)$ corrections coming from the 1-loop QCD vertex correction Cline:2013gha , which yields the following correction factor Drees:1990dq to the fermionic expressions (7)

[TABLE]

and the factor of three comes from the number of color charges. For better accuracy, we use the higher order perturbative results from the HDecay code Djouadi:1997yw for $m_{S}>2.5{~{}\rm GeV}$.

Metastable mesons, such as $\pi^{\pm}$ and $K^{\pm},\bar{K}^{0},K^{0}$ are “important” decay products, as they can participate in the charge-exchange reactions with nucleons and shift the $n-p$ balance, hence affecting the whole nucleosynthetic chain. In the mass range where the perturbative QCD calculations are no longer valid, we base our baseline calculations on Ref. Bezrukov:2009yw . The scalar-pion interaction can be extracted from the low-energy expansion of the trace of the QCD energy-momentum tensor (see for e.g.. Voloshin:1985tc ; Leutwyler:1989xj ) by integrating out the three heavy quarks and using chiral perturbation theory on the remainder, yielding the effective Lagrangian Bezrukov:2009yw

[TABLE]

where we have inserted the SM numerical values for the number of heavy quarks and the first coefficient of the QCD beta function. This leads to decay width to charged pions

[TABLE]

This result is however not applicable far above the pion threshold, as final-state resonances would drastically affect this prediction. Instead, we use the pion and kaon decay width described in Ref. Donoghue:1990xh , where the authors matched the next-to-leading order corrections of the low-energy theorems to the dispersion results from the $\pi\pi$ phase-shift analysis above 600 MeV from the CERN-Munich group Hyams:1973zf . The photon decay channel is added with the prescription detailed in Ref. Spira:1997dg . Finally, there is a gap for $1.4~{}{~{}\rm GeV}<m_{S}<2.5~{}{~{}\rm GeV}$ where no analytical treatment is entirely trustworthy, as this includes new resonances strongly coupled to $\eta\eta$ and other potential hadronic cannels. We simply follow Ref. Bezrukov:2009yw and interpolate between the two regimes, under the assumption that there is no order of magnitude deviation in this mass range. The branching ratios and the lifetime for $\theta=10^{-6}$ are displayed in Fig. 1.

As an alternative decay spectrum model, we also display the perturbative spectator approach Gunion:1989we ; McKeen:2008gd ; Alekhin:2015byh , where the relative decay width above the kaon threshold are given by

[TABLE]

with $\beta_{i}=\sqrt{1-4m_{i}^{2}/m_{S}^{2}}\Theta(m_{S}-2m_{i})$, $\Theta$ being the step-function, and we adopt the running of $s$ quark mass following Ref. Spira:1997dg . The pion contribution is kept as in equation (10) and then we use the HDecay output at the $c$-quark threshold and above to match our baseline model.

For $m_{S}$ of several GeV and heavier, decays with final state nucleon-antinucleon pairs are possible. Even though the branching to such states are generally lower than 10%, the effect on BBN can be quite significant, and therefore these are by far the most important channels for $\tau_{S}\gtrsim 1\;\sec$. On top of direct and for the most part subdominant contributions from $S\to\bar{n}n,...$, we need to take into account the (anti-)nucleon states that emerge from the hadronization of the quark decay products and heavy $B$-meson fragmentations.

II.2 Cosmological metastable abundance

After the temperature drops below $m_{S}$, the interaction of $SS$ pairs with the SM shifts towards the annihilation, resulting in an intermediate (metastable) population of $S$ bosons. In the mass range that we consider, the $S$ annihilation is dominated by the $s$-channel reactions $SS\to h^{*}\to XX$, where on the receiving end are the pairs of the SM states $XX$ created by a Higgs-mediation process. The annihilation cross section $\sigma v$ generically takes the form

[TABLE]

This formula recast the rate in terms of a Higgs width $\Gamma_{{\rm SM}}^{m_{h}\to\sqrt{s}}$ with a fictitious mass of $\sqrt{s}$. This form encompasses both perturbative and non-perturbative channels in the $h^{*}$ decay rate (with the substitution $m_{h}^{*}\to\sqrt{s}$), which we have described above. In the standard WIMP freeze out paradigm, a DM particle freezes out at $T_{\rm f.o.}\sim m_{DM}/20$, $\langle\sigma v\rangle$ is simply the nonrelativistic limit $\sigma v(\sqrt{s}=2m_{DM})$ and the relic density can be conveniently approximated as $\Omega_{DM}h^{2}\sim 0.11\times 1{\rm pb}/\langle\sigma v\rangle$. This result emerges as a solution to the Boltzmann equation111We use the standard variable definitions, where $Y=n_{S}/s$ is the $S$ abundance normalized on the entropy density $s$, $x=m/T$ is the dimensionless inverse temperature, $H$ is the Hubble rate, $h_{\rm eff}$ is number of entropic relativistic degrees of freedom and $Y_{eq}$ is the normalized thermal equilibrium $S$ number density. Steigman:2012nb

[TABLE]

when the freeze out occurs in the exponentially falling region of the equilibrium density $Y_{\rm eq}(T)$. For a much smaller annihilation cross section, $\langle\sigma v\rangle\ll 1\;{\rm pb}$, $Y$ departs from the equilibrium value earlier, possibly near the relativistic plateau $Y_{\rm eq}=n_{\rm eq}/s\to 45\zeta(3)/2\pi^{4}h_{\rm eff}(T)$ for $x\ll 1$. Since the nonrelativistic annihilation cross section in the minimal Higgs portal model ranges from $10^{-3}$ to $10^{-14}$ pb for $m_{S}\sim 1{~{}\rm MeV}-60$ GeV and $Br(h\to SS)\sim 0.1-0.001$, we numerically integrate equation (13) to determine the metastable $S$ abundance. The results are shown in Fig. 2, normalized to the baryon number density for a more intuitive interpretation of its impact on BBN in the following section.

For $m_{S}\simeq m_{h}/2$, the $\sigma v$ cross section evaluated at $s=4m_{S}^{2}$ is a poor approximation, as it fails to capture the strong energy dependence of the cross section near the resonance at $\sqrt{s}=m_{h}/2$ Griest:1990kh . The sharp drop in the abundance above $m_{S}\sim 45{~{}\rm GeV}$ is due to the resonant contribution to the thermally averaged cross section, leading to a delayed freeze out and drastic decrease in metastable $S$ abundance. Our numerical results agree with the semi-analytic treatment of Ref. Cline:2013gha . For very light $m_{S}$, one can see that the freeze out abundances are large, and the relative spread between different input values of $Br(h\to SS)$ gets smaller, as the annihilation cross section becomes very small and the freeze out happens in the semi-relativistic regime $x_{\rm f.o.}\sim\mathcal{O}(1)$ and asymptote to the $Y_{\rm eq}$ relativistic plateau for small $m_{S}$. The only difference at the lightest masses is from $Y_{\rm eq}^{rel}\propto 1/h_{\rm eff}(T)$. Since $h_{\rm eff}$ is a monotonic function of temperature, weaker annihilation cross sections freeze out earlier, at a higher temperature, thus yielding smaller abundances (as seen in the $m_{S}=5{~{}\rm MeV}$ curves in Fig. 2). This is in contrast with the standard freeze out in the non-relativistic regime, with final abundances inversely proportional to the cross section. We note in passing that the strong-interaction-related uncertainty “propagates” outside the $m_{S}\sim 2m_{\pi}-2m_{c}$ window. For example, because of the relativistic freeze out, for $m_{S}$ smaller $2m_{\pi}$ the hadronic channels may turn out to be important.

III Big Bang Nucleosynthesis

The formation of light nuclei is one of the earliest probes of NP in cosmology along with far less certain constraints imposed by the inflationary framework. BBN is well-understood within SM physics, and its outcome agrees with observational data for 4He and D. 7Li has an outstanding factor of $\sim 2-3$ discrepancy between theory and observations Cyburt:2015mya , with the caveat that the observed abundances may have been affected by stellar evolution. Nevertheless, the overall success over a wide range of abundances can be used to constrain various types of NP Pospelov:2010hj .

The initial BBN stage is the neutron-proton ratio $n/p$ freeze out. Maintained in equilibrium by electroweak interactions at high temperatures, the neutron abundance follows $n/p\sim e^{-Q/T}$, where $Q=m_{n}-m_{p}-m_{e}\simeq 1.293{~{}\rm MeV}$, until the epoch when the weak processes decouple around temperatures of 0.7 MeV. The outcome, $n/p\simeq 1/6$, is quasi-stable, decreasing to $n/p\simeq 1/7$ at the end of the “deuterium bottleneck”. The latter terminology is used to indicate a much delayed onset of nuclear reactions controlled by a relatively shallow $n-p$ binding energy. Once the Universe runs out of photons that can efficiently dissociate deuterium, the bulk of the nucleosynthetic reactions occurs at $t_{\rm deut}\sim 200$ seconds. 4He has a large binding energy per nucleon, and the reactions leading to it are less Coulomb-suppressed than for heavier elements. Consequently, most neutrons end up in the final 4He abundance (expressed in mass fraction from the total baryon mass) $Y_{p}\simeq 2(n/p)/\left(1+n/p\right)\simeq 0.25$.

Traces of neutrons and incomplete nuclear burning of $A=2,~{}3$ nuclei light nuclei result in the left-over abundances of 3He and D. Beyond the 4He atomic number, the deepest bound nucleus is 12C, but its formation is completely suppressed since it would need to be produced by a triple 4He collision. The $2\to 2$ reactions $p+{}^{4}$He and 4He + 4He are also ineffective at producing heavier nuclei as the $A=5$ and $A=8$ elements are all unstable. The only remaining possibilities are ${}^{4}{\rm He}+{}^{3}{\rm He}\to{}^{7}{\rm Be}+\gamma$ followed by a $\beta$ decay to yield ${}^{7}{\rm Li/H}\sim\mathcal{O}(10^{-10})$ and 6Li formed at the ${}^{6}{\rm Li/H}\sim\mathcal{O}(10^{-14})$ via 4He-D fusion. For the problem at hand - the determination of the upper limit on the $S$ lifetime - few of these details matter. This is because of relatively large metastable abundances affecting the earliest stages of nucleosynthesis, primarily via the $n/p$ ratio.

III.1 Neutron Enrichment

Ample abundances of $S$ particles ($n_{S}\sim 10^{2}-10^{9}\times n_{b}$) flood the Universe with final state mesons and nucleons that in turn could spoil the final light nuclei abundances. For example, at temperatures $T\sim 0.5$ MeV, the protons are $\sim 6$ times more abundant than neutrons, but this ratio can be easily changed due to meson-induced charge exchange reactions. At these temperatures, the probability of $p\to n$ conversion from charged pions is

[TABLE]

It is then clear that injection of $O(10^{3})$ mesons per nucleon at these temperatures can drastically increase the $n/p$ freeze out abundance. Similarly, direct baryonic injection of $n\bar{n}$ and $p\bar{p}$ will have a similar effect on the $n/p$ ratio. On the other hand, if $S$ decays happen before the $n/p$ freeze out, the additional $p\to n$ conversions would not be as efficient, being washed out by the on-going weak interaction conversions.

The limit of the exclusion region in the $Y_{S}/\tau_{S}$ parameter space ($Y_{S}\equiv n_{S}/n_{b}$ from now on) is determined by solving the Boltzmann equation with the injection of charge exchange inducing particles. Given that the abundances of $S$ particles are large, the main constraints can be derived from the $n/p$ freeze out ratio. To that effect, we would not need a complete BBN framework, but only a subset of the whole code that deals with $n\leftrightarrow p$ conversions. We follow the semi-analytic treatment by Mukhanov Mukhanov:2003xs , that approximates $n\leftrightarrow p$ weak conversion rates by a few integrals over thermal distributions, and assumes a “step-like” disappearance of charged leptons below $T=m_{e}$,

[TABLE]

where $g_{a}\simeq 1.27$ is the standard nucleon axial-vector coupling, $Q=m_{n}-m_{p}-m_{e}\simeq 1.293{~{}\rm MeV}$, and $G_{F}$ is the Fermi constant. The reverse reaction rates are found by detailed balance. We evaluate $J$ numerically and solve for the electron-neutrino temperature $T_{\nu}$ by entropy conservation, assuming a $\nu_{e}$ decoupling temperature of 2 MeV, which reproduces the correct entropy degrees of freedom at lower temperature Srednicki:1988ce . It is then straightforward to solve numerically the differential equation for $X_{n}=n_{n}/n_{b}$,

[TABLE]

where the last term represents the neutron decay with $\Gamma_{n}^{-1}=880~{}{\rm sec}$. This equation is approximately valid until the rapid switch-on of the nuclear reaction rates at the end of the deuterium bottleneck. Within this approximation, one can determine the final temperature where the equation is valid by starting with $X_{n}=1/2$ at early times, and solving for the deuterium bottleneck temperature by imposing $Y_{p}=2X_{n}(T_{\rm deut})=0.25$. This results in $T_{\rm deut}\simeq 0.068{~{}\rm MeV}$ or $t_{\rm deut}\simeq 276$ sec. We take this approximation as our baseline SBBN model, which is then modified by the inclusion of extra sources and sinks for $n,~{}p$, and new $n\leftrightarrow p$ reactions. To constrain the parameter space of a species decaying into charged mesons or baryons, we proceed by solving the Boltzmann equation that includes new interactions. We will require that $Y_{p}$ does not deviate from SBBN by more than 4%,

[TABLE]

which is a rather generous allowance for the errors, considering the tight observational constraints on primordial helium abundance Cyburt:2015mya . Consequently, it will result in conservative limits of $\tau_{S}$.

III.1.1 Meson-Mediated mechanism

Only long-lived mesons have an opportunity to interact with the baryon bath and induce proton-neutron conversions. As such, only $\pi^{\pm}$, $K^{\pm}$ and $K_{L}$ have lifetimes in excess of $\tau\sim 10^{-8}$ seconds, and can induce $p\leftrightarrow n$ via strong interactions. For temperatures relevant for the $n/p$ freeze out, the density of charged leptons is very high, and mesons are efficiently stopped by the primordial plasma. We assume that they are efficiently thermalized, and take the relevant pion-induced reactions at threshold Reno:1987qw ; Pospelov:2010cw ($c=1$),

[TABLE]

The reverse reactions are irrelevant due to the short lifetime of $\pi^{0}$’s and the need for non-thermal $\gamma$’s of $\sim 140$ MeV energy. The $\pi^{-}$ reactions are to be added to the r.h.s of Boltzmann equation (18) via the additional term

[TABLE]

and similarly for the $\pi^{+}$ reactions. The ambient population of injected pions from a $S$ decay with $Br(S\to\pi^{+}\pi^{-})=\xi_{\pi^{\pm}}$ is $n_{\pi^{\pm}}^{\rm inj}\simeq\xi_{\pi^{\pm}}\Gamma_{S}\tau_{\pi^{\pm}}Y_{S}n_{b}(T)e^{-t\Gamma_{S}}$, $t\simeq 2.42\sec\;({\rm MeV}/T)^{2}/\sqrt{g_{\star}}$ and the thermal cross section are taken at their threshold value $\langle\sigma v\rangle_{np}^{\pi^{+}}=(\sigma v)_{np}^{\pi^{+}}$. Reactions with pairs of charged particles in the initial states, such as $\pi^{-}p$, will be somewhat enhanced due to the Coulomb attraction, which provides a small but non-negligible correction. We account for it following Ref. Pospelov:2010cw .

The implementation of the charged kaons reactions is similar to the pion case, but the dominant reactions are rather different. The direct charge exchange between neutral and charged kaons is

[TABLE]

with similar cross section for a charge-conjugated reaction, $K^{0}p\to K^{+}n$. For neutral kaons, the effects induced by $K_{L}$ are the most important, and we use $\sigma(K_{L}n\to K^{-}p)\simeq\frac{1}{2}\sigma(\bar{K}^{0}n\to K^{-}p)$ and (25) to find $(\sigma v)^{\bar{K}^{0}}_{pn(K^{-})}\simeq 4.5{~{}\rm mb}$. Additionally, efficient reactions can also proceed via $s$-quark being incorporated inside a hyperon that subsequently decays into $p/n$ + $X$. The inclusive threshold cross section found by weighting each hyperon with their branching ratios to $p/n$ are Pospelov:2010cw

[TABLE]

Notice the absence of corresponding hyperon reactions initiated by $K^{+}$ due to the presence of anti-$s$ quark.

Representative examples of $X_{n}(T)$ evolutions and the sensitivity to $\xi_{\pi^{\pm}}Y_{S}/\tau_{S}$ parameter space are shown in Fig. 3. Left panel displays significant modifications to the evolution of neutron abundance at $\tau_{S}=$0.05 seconds with adjustable initial abundance, yielding $\Delta Y_{p}=0.01$. The departure from $X_{n}=0.5$ at high temperatures is clearly visible. (In fact, for short $\tau_{S}$, the kaon injection channel at early times leads to a shift of the equilibrium value of $X_{n}$ to $(\sigma v)^{K^{-}+K_{L}}_{pn}/((\sigma v)^{K^{-}+K_{L}}_{pn}+(\sigma v)^{K^{-}+K_{L}}_{np})\simeq 0.45$.) As the temperature lowers, the Coulomb-enhanced reaction becomes stronger. For meson injection, these reactions enhance the $p\to n$ conversion, keeping $X_{n}$ away from the SBBN value. Right panel gives a boundary of the exclusion regions for different injection modes. In addition to the already described channels, charged kaons also give rise to a population of secondary charged pions that should also be included in the analysis of $p\leftrightarrow n$ transitions. Since the constraints are already stronger than for the charged pion case, we neglect this effect, which leads to more conservative bounds.

III.1.2 Direct baryonic injection mechanism

If $S$ is heavy enough, the end-products after hadronization of the primary decay products (e.g. $b$ or $c$ quarks) may contain baryons. Since $S$ has no baryon number, one should expect an equal number of baryons and anti-baryons in the final states. Therefore, one should expect the injection of $n\bar{n}$, $p\bar{p}$, $\bar{n}p$ and $p\bar{n}$ pairs, as well as (in principle) baryonic states with higher multiplicities. The hadronization process and decay of heavy quarks produce much more light mesons than baryons, and a complete analysis must include a Monte Carlo study of the hadronization process (see Ref. Kohri:2001jx for benchmarks of heavy unstable particles decaying into 2 hadronic + 1 leptonic jets in the early BBN epoch). Assuming that the heavy quarks inside baryons decay due to the “main” weak decay sequence, $b\to c\to s\to u$, one should also expect a somewhat large number of the final states with a proton or anti-proton over neutron or anti-neutron. We will tune the branching models of $S$ to available particle data on proton production, and take $N_{n}=\kappa N_{p}$ and $N_{\bar{n}}=\kappa N_{\bar{p}}$. Furthermore, due to a more frequent appearance of up-quark at the end of the decay chain, we would take $\kappa\simeq 0.5$ on average.

As in the case of mesons, the thermalization of baryonic decay products is quick (see e.g. Kawasaki:2004qu ). As a baryonic pair is created in the decay, the baryon is added to the existing population of $n$ or $p$. The anti-baryon will, however, annihilate with either $p$ or $n$ and dissipate into lighter mesons. If it annihilates with its own antiparticle, there is no net change in $n/p$, but an annihilation with the other species induces a net $n-p$ change. The probability $P_{k\bar{l}}^{i\to j}$ of a net charge exchange $i\to j$ from a $k\bar{l}$ injection is simply given by the weighted annihilation rates

[TABLE]

and similarly for the $n\bar{p}$ and $n\bar{n}$ injections. The baryonic annihilation rates are given by  Reno:1987qw

[TABLE]

where the $p\bar{p}$ has the low-$v$ Coulomb correction $C(v)$. The implementation of these processes in the Boltzmann equation then require additional terms

[TABLE]

As before, the outcome is displayed in Fig. 3. Again, for short $S$ lifetimes and large $Y_{S}$, the large numbers of injected particles completely dictates the early $X_{n}$ value. The constraint on $Y_{S}$ goes up more sharply in the short $S$ lifetime limit. There is a significant dependence on $\kappa$ for $\tau_{S}\gtrsim 0.1$ sec, which is washed out by the SM electroweak interactions at earlier times. If we take the extreme limit $\kappa\to 0$, no neutrons are injected and the $p\bar{p}$ pair can only further decrease the $n/p$ ratio, thus constrained by the lower $Y_{p}$ limit 0.24. On the other hand, a symmetric injection $\kappa=1$ enhances the $n/p$ ratio as the anti-baryon mostly annihilates on protons, more abundant than neutrons by a factor of $\sim 6-7$ after the standard $n/p$ freeze out. For $\kappa\gtrsim 0$, the final $Y_{p}$ can either be increased or decreased, depending if the $S$ particles decay away before or after the displaced $X_{n}$ equilibrium crosses the SBBN $n/p$ freeze out curve. As shown for $\kappa=0.2$ in Fig. 3, there is a $Y_{p}=0.24$ exclusion island at low lifetimes and larger lifetimes are constrained by $Y_{p}=0.26$. We use $\kappa=0.5$ as our baryonic injection constraint benchmark.

III.1.3 Muon-Mediated mechanism

Muon injection physics differs from the previous scenarios of meson and baryon injection. The direct charge-exchange is through the weak force, as opposed to the strong force in the other cases, and is completely negligible over the lifetime of the muon. Instead, the reactions can happen via the energetic neutrinos emitted by the muon decays. The case for muon injection after $t\sim 100\;\sec$ has been covered in Ref. Pospelov:2010cw , to which we refer the reader for details. Assuming stopped muons, the authors solved for the injected neutrino energy spectrum, including redshifting and averaged over flavour oscillations, to be integrated in the $n-p$ conversion rate. At earlier times, we know background neutrinos are coupled to $e^{\pm}$ down to $T\simeq 2{~{}\rm MeV}$, and energetic injected neutrinos must accordingly deplete their energy efficiently as well. Summing over the possible interactions with the background neutrinos and $e^{\pm}$ Dolgov:1997mb , the collision rate of an injected electron-neutrino with the bath is given by

[TABLE]

$g_{\rm L}=1/2+\sin^{2}\theta_{\rm w}$, $g_{\rm R}=\sin^{2}\theta_{\rm w}$, while $\eta(T)=1$ for $T\gtrsim m_{e}$ and exponentially falls to 0 at lower temperatures. We follow the implementation of Ref. Pospelov:2010cw and correct for the removal of energetic neutrinos by adding an effective collision lifetime in the neutrino energy distribution (normalized on $n_{b}$)

[TABLE]

where $F_{e}$ is the distribution at injection time $T_{1}$, averaged over flavour oscillations. The charge-exchange rate to be inserted in the Boltzmann equation (18) is

[TABLE]

and similarly for the reverse $np$ direction. The resulting constraints are shown in Fig. 3. Our results lean on the conservative side on a few assumptions. For simplicity, we assumed one collision for the neutrino thermalization, instead of following energy degradation over a shower of multiple interactions. Moreover, we took the collision time of the electron-neutrino, even though there are muon-neutrino states in the oscillations. Since $\Gamma_{\rm coll}^{\nu_{e}}>\Gamma_{\rm coll}^{\nu_{\mu}}$, we over-estimate the actual collision time and the overall conversion rate should be slightly larger.

III.1.4 Meson injection from residual annihilations

In addition to its decay products, $S$ can also inject particles in the cosmic medium via $SS$ annihilations to charged pions. The injected pions interact with the cosmic medium in the same way as from $S$ decays described above. The Boltzmann equation takes the addition term (24), with the injected pion density now given by

[TABLE]

where $\langle\sigma v\rangle_{\pi^{+}\pi^{-}}$ is the non-relativistic annihilation cross section $\sigma v(2m_{S})$ as per Eq. (12), rescaled by the pionic branching ratio at $\sqrt{s}=2m_{S}$. The $n_{\pi^{\pm}}^{\rm ann}\propto n_{S}^{2}\propto T^{6}$ dependence imply a much stronger impact at high energies, enforcing the displaced initial condition $X_{n}^{i}\simeq 0.47$. As $S$ decays away, its impact on $X_{n}(T)$ is even more rapidly exponentially suppressed and its constraints are less stringent than decays at very short lifetimes. The bounds from annihilation are given in the $Y_{S}^{2}\langle\sigma v\rangle_{\pi^{+}\pi^{-}}-\tau_{S}$ parameter space and displayed in Fig. 4.

III.2 Energy density requirements

The resultant BBN abundances depend on the nuclear reaction rates and how efficient they are as the Universe expands. One by one, the reaction rates drop out of equilibrium, as the Universe expands and cools. If the Hubble rate is increased due to a large energy density locked in a dark sector, the active reaction time would shorten, potentially spoiling the SBBN results. For our study, the most important effect is the change of the Hubble rate during the $n/p$ freeze out, which again affects $Y_{p}$. However, we can also use as a constraint a well-measured quantity in cosmology, the total energy density carried by neutrinos.

The neutrinos decouple from thermal processes at $T\sim 2{~{}\rm MeV}$. If the decaying particle is heavy and does not decay into neutrinos, it will reheat electron-photon fluid with respect to the neutrinos, decrease $T_{\nu}/T_{\gamma}$ and equivalently lower $N_{\rm eff}$. The Planck collaboration measured $N_{\rm eff}=3.04\pm 0.33$ at $2\sigma$, including their CMB results and external cosmological data Ade:2015xua , which imposes $N_{\rm eff}>2.71$ as a lower bound.

The energy densities and Hubble rate form a closed system of differential equations

[TABLE]

where we have assumed a non-relativistic $S$ and omitted the variation in relativistic degrees of freedom. Assuming step-like decoupling and changes in relativistic degrees of freedom, the $T$ evolution separates into 3 regions. For $T>T_{\nu}^{\rm decoup}$, neutrinos are in equilibrium with the electromagnetic bath and $\rho_{S}$ is injected equally in $e^{\pm}$’s, $\nu$’s and $\gamma$’s. For $T_{\nu}^{\rm decoup}>T>T_{m_{e}}$, the neutrinos are simply redshifted while the electron-photon bath is heated by the $S$ decays. For $T_{m_{e}}>T$, electrons become non-relativistic and transfer their entropy to photons, additionally heating the photon bath compared to the neutrino bath.

If $S$ does not dominate the energy density of the Universe before its decay, we can write $\rho_{S}=\delta_{S}\rho_{rad}^{SM}$, $\rho_{rad}=\rho_{rad}^{SM}(1+\delta_{rad})$ and expand (37) around the $\delta$ perturbations to solve the system analytically. At linear order, we find the solutions

[TABLE]

where $\Gamma_{3/2}$, $\Gamma_{5/2}$ are incomplete Gamma functions and the integration constants $c_{S}$, $c_{rad}$ are set to have $\rho_{S}=m_{S}n_{S}$ and $\rho_{rad}=\rho_{rad}^{SM}$ at some early time $\Gamma_{S}t\ll 1$. After the neutrinos decouple, the injected energy is distributed to the photon-electron bath and its energy density departs for the neutrino bath

[TABLE]

where $\tilde{g}_{i}\equiv g_{i}/(g_{\gamma+e}+g_{\nu})$ is the fraction of relativistic degrees of freedom of each bath, $t_{\nu}^{\rm decoup}$ the neutrino decoupling time and

[TABLE]

Finally, after the electrons become non-relativistic, they effectively transfer their entropy to the photon bath. Assuming an instantaneous transition, entropy continuity implies an increase of energy density by a factor of $\delta=\left(g_{e}+g_{\gamma}\right)^{1/3}/g_{\gamma}^{1/3}=(11/4)^{1/3}$. Matching boundary conditions, the energy densities at late times are

[TABLE]

with $\alpha=1/(\delta\tilde{g}_{\gamma+e}+\tilde{g}_{\nu})$ and $c$ a boundary condition that is irrelevant in the $t\to\infty$ limit.

The temperature-time dependence is found via $\rho_{rad}(t)=\pi^{2}g_{\star}T^{4}/30$. Since the neutrino interaction rate scales as $\Gamma_{\nu_{e}}\sim T^{5}$, we find the neutrino decoupling time in the modified cosmology by equating $(T_{\nu}^{\rm decoupl})^{5}/H(T_{\nu}^{\rm decoupl})=(T_{\nu}^{0})^{5}/H_{0}(T_{\nu}^{0})$, with $H$ the perturbed Hubble rate and $T_{\nu}^{0}$ the neutrino decoupling temperature in the SM. In the Maxwell-Boltzmann approximation, $T_{\nu}^{0}=2{~{}\rm MeV}$, but thermal refinements in the interaction rates and phase space tend to yield a lower value $T_{\nu}^{0}=1.4{~{}\rm MeV}$ Dolgov:2002wy .

Then, at $\Gamma_{S}t\gg 1$, we can evaluate $T_{\nu}/T_{\gamma}$ and find

[TABLE]

to constrain energy injection into electrons. We display in Fig. 5 the departure from $N_{\rm eff}=3$ as a function of time for $\tau_{S}=0.1\sec$ and the two neutrino decoupling temperature benchmarks. The limits are also shown in units of stored energy density $m_{S}n_{S}/n_{b}$, where $m_{S}$ is in MeV. If the $S$ decay happens after the neutrino decoupling, all energy is deposited in the photon bath and the result is independent of our choice of $T_{\nu}^{0}$. If decays happen earlier, the photon and neutrinos are potentially still coupled and the energy emitted in $S$ decays only influences $N_{\rm eff}$ after decoupling. As such, the constraints has a $t\propto(1/T_{\nu}^{0})^{2}$ dependence. We adopt the conservative side, $T_{\nu}^{0}=1.4{~{}\rm MeV}$, as our bounds on the $m_{S}-\tau_{S}$ parameter space. Notice the constraints for $SS$ annihilations to pions are much stronger and will be dominant when the pionic annihilation channel is open, i.e. for $m_{\pi}<m_{S}<2m_{\pi}$.

III.2.1 Energy injection partitioned between photon and neutrino baths (e.g. muon injection)

The case for muon injection is somewhat interesting as its decay products, neutrinos and electrons, clearly thermalize in the two different baths, once everything is decoupled. Both $T_{\gamma}$ and $T_{\nu}$ will rise, but since the 2 neutrinos carry more energy than the electron for a muon decay, we expect a rise in $N_{\rm eff}$. More precisely, we solve a similar set of equations as (37), except the photon bath absorbs a $\xi$ proportion of the $S$ decay energy and the neutrino bath gets the remaining $(1-\xi)$ portion. Before neutrino decoupling, the radiation bath evolves as in equation (38). Each decay product carries on average the energy Pospelov:2010cw

[TABLE]

After neutrino decoupling, the energetic neutrinos can still collide with the ambient electrons until $\Gamma_{{\rm coll}-e}^{\nu_{e}}<H$, where $\Gamma_{{\rm coll}-e}^{\nu_{e}}$ is the collision rate with electrons only, the $T_{\gamma}$-dependent term in equation (III.1.3). Then, the energy distributed to the photon bath separates into two regimes

[TABLE]

where the muon-neutrino collision term is given by

[TABLE]

Following the same procedure as before, we find

[TABLE]

with $t_{\rm coll}$ found by solving $\Gamma_{{\rm coll}-e}^{\nu_{e}}=H$. The physics is constrained by $N_{\rm eff}<3.37$. The time-dependence of the departure from $N_{\rm eff}=3$ is shown for $\tau_{S}=0.2\;\sec$ and the two choices of $T_{\nu}^{0}$ in Fig. 5. The corresponding constraints on the maximal stored energy for a given lifetime are shown on the right. For comparison with the muon-induced $Y_{p}$ bound, we display the curve for $m_{S}=250{~{}\rm MeV}$ from neutron enrichment in purple. Independently from the choice of $T_{\nu}^{0}$, the bounds from $Y_{p}$ $\mu$-injection are most constraining for $\tau_{S}\gtrsim 0.2\;\sec$ while the annihilation to $\pi^{+}\pi^{-}$ provides the dominant constraint in the entire $2m_{\mu}<m_{S}<2m_{\pi}$ range.

III.3 Late-time energy injection

In the example of the $S$ particles coupled through the Higgs portal, the most stringent constraints on lifetime come from the considerations of $n/p$ freeze out. In other models, with additional channels of annihilation that can suppress metastable abundances, the constraints on lifetime would not be as stringent, and would mostly come from the considerations of late energy injection. For completeness, we also discuss these constraints here. Modification of BBN by unstable particles with lifetimes in excess of 200 seconds has been considered in detail, both through hadronic Reno:1987qw ; Kawasaki:2004yh , electromagnetic Cyburt:2002uv or combined Kawasaki:2004qu ; Jedamzik:2006xz ; Pospelov:2010cw energy cascades.

Hadronic injection after $t\gtrsim 200$ seconds is most efficient at modifying the final yields of the less abundant light nuclei D, 3He, 6Li, and 7Li. After most of 4He has been synthesized, the BBN enters the regime ($T\sim 50$ keV) when neutrons are rare, $O(10^{-5})$ or so, yet their abundances are critical in determining the final abundance of deuterium. At that stage, any additional neutrons brought into the system through external processes such as heavy particle decays lead to the increase of the deuterium abundance. (Incidentally, it also leads to the suppression of 7Be and consequently of 7Li Reno:1987qw .) The increase of D production can be exacerbated by the hadro-dissociation of 4He in the process of slowing down of injected hadrons. Additional production of 3He through spallation can also affect the 3He/D ratio Ellis:2005ii . Secondary and tertiary processes may also generate 6Li and 9Be Dimopoulos:1987fz ; Pospelov:2010kq . Detailed studies of the ensuing constraints Jedamzik:2006xz show strong sensitivity to hadronic (mostly nucleonic) decays of metastable particles with lifetimes in the hundreds of seconds and longer, and initial abundances comparable or even smaller than that of baryons. In recent years, these constraints have only got stronger, primarily due to steady observational progress in determination of primordial D/H Cooke:2013cba .

If for some reason, hadrons and specifically nucleons are absent from the decay chains, the abundances of light elements can be modified by the late injection of electromagnetic energy. At early times this mechanism is inefficient, as radiation quanta with energy in excess of nuclear binding are quickly energetically degraded by ambient plasma. The photo-dissociation therefore sets in at late times leading to a suppression D ($t\gtrsim 10^{4}$ seconds) and additional production of 3He for $t\gtrsim 10^{6}$ seconds. Since typically 45% of hadronic energy injection is dissipated electromagnetically in the hadronization cascade Jedamzik:2006xz , the late-time energy injection constraints on a heavy particle are dominated by the electromagnetic reactions in the BBN network.

IV Results

We are now in a position to perform a scan in parameter space of the minimal Higgs model, constrained by the consistency with BBN. In Fig. 6, we display the parameter space, both in the lifetime and an effective decay length $L_{\rm dec}=c\tau_{S}\beta_{S}(E_{S}/m_{S})$. We assume an average $E_{S}$ of 200 GeV, from a Higgs typically boosted at 400 GeV at the LHC. The resulting constraints, along with the assumptions considered in each mass range are described below.

• •

Region A $2m_{e}<m_{S}<2m_{\mu}$ : The constraint comes from the decrease in $N_{\rm eff}$ with the entropy dump in the SM bath after neutrino decoupling. We take the neutrino decoupling temperature to be $T_{\nu}^{0}=1.4{~{}\rm MeV}$ as a conservative limit.

• •

Region B $m_{\pi}<m_{S}<2m_{\pi}$ : This region is dominated by the $SS$ annihilation to $\pi^{+}\pi^{-}$. We also derived the same constraint as region A from $N_{\rm eff}$ up to $m_{S}=2m_{\mu}$, in addition to the raised $N_{\rm eff}$ from decays into muons in the $2m_{\mu}<m_{S}<2m_{\pi}$ and the $Y_{p}$ constraints from $S$ decaying into muons. They all yield weaker bounds, of $\tau_{S}>0.3\sec$ or longer.

• •

Region C $2m_{\pi}<m_{S}<2m_{K}$ : The abundance $Y_{S}$ weighted by the pion branching ratio constrains the region via direct charged pion decays. We assume $2/3$ go into charged pions and $1/3$ is radiated away in $\pi^{0}$.

• •

Region D $2m_{K}<m_{S}<1.4{~{}\rm GeV}$ : The abundance $Y_{S}$ weighted by the kaon branching ratio constrains the region via direct charged kaons decays. We assume $1/2$ go into charged kaons and $1/2$ into $K^{0}\bar{K}^{0}$. Only half of the neutral kaons survive as $K_{L}$, creating similar in numbers metastable populations of $K_{L}$, $K^{+}$ and $K^{-}$.

• •

Region E $1.4{~{}\rm GeV}<m_{S}<2m_{D}$ : By strangeness conservation, we assume that all $s$-quarks yield a kaon, half charged and half neutral. Since we do not have model-independent branching ratios of $S$ in this mass regime, we vary the description according to the assumptions in each decay model. For the baseline model, we assume that 100% decays to the kaons and apply our kaon injection constraints. For the perturbative spectator model, the kaon branching ratio is given by (11), with non-negligible contributions from decays to pions, muons and eta mesons, resulting in weaker bounds until the $c$-quark threshold. At $m_{S}=m_{c}$ the hadronic modelling dependence largely goes away.

• •

Region F $2m_{D}<m_{S}<2m_{b}$ : We utilize the branching fractions of $c\bar{c}$ from $e^{+}e^{-}$ at $\sqrt{s}=10.5{~{}\rm GeV}$ into $D$-mesons from Ref. Lisovyi:2015uqa and weight each channel by its inclusive $K^{\pm}$ branching ratios to find a hadronization yield of 0.63 $K^{+}K^{-}$ pair per $S$ decay into $c$-quarks. Rescaled by $Br(S\to c\bar{c})$, same constraints from kaon injection apply. Above the $2m_{\Lambda_{c}}$ threshold, a $c\bar{c}$ typically forms a $c$-baryon with a 0.06 probability Lisovyi:2015uqa , which then hadronizes to $p$ or $n$. We find this constraint weaker than the kaons injection and use the $K^{+}K^{-}$ result across this entire range.

• •

Region G $m_{S}>2m_{b}$ : The main decay channel here are pairs of $b\bar{b}$ quarks. The charged pion, charged kaons and proton multiplicities in the $b\bar{b}$ decay of a $Z$ boson are measured to be $18.44\pm 0.63$, $2.63\pm 0.14$ and $1.00\pm 0.08$ respectively by the ALEPH collaboration Barate:1997ty . We assume the ratio holds in the hadronization of lower centre-of-mass decays into $b\bar{b}$ and scale by the mean charge multiplicity fit Sarkisyan:2015gca

[TABLE]

where $s_{0}=1{~{}\rm GeV}^{2}$. This fit agrees well in both $e^{+}e^{-}$ and $p\bar{p}$ collisions between $\sqrt{s}\sim 2{~{}\rm GeV}-2{~{}\rm TeV}$. This gives us an estimate for the baryon injection of the $b\bar{b}$ branching fraction of $S$. We further assume 50% smaller injection of $n(\bar{n})$ to utilize our baryon injection constraints. The accompanying pions and kaons also independently yield comparable constraints, not shown in the figure.

V Discussion

We have considered, in some detail, constraints on the lifetimes of the scalar particles, coupled to the Higgs portal via a minimal set of couplings. To stay relevant for the LHC, we have concentrated on $m_{S}<m_{h}/2$ case, that allows pair-production of $S$ states in the decay of Higgs bosons. The same coupling is responsible for the cosmological depletion of $S$ particles, leading to their metastable abundance in the early Universe.

We find that throughout almost the whole mass range considered in this work, $2m_{\mu}<m_{S}<m_{h}/2$, the constraints on the lifetime of $S$ particles are stronger than $0.1$ seconds. Moreover, the results have a relatively mild dependence on the $Br(h\to SS)$. The reason for that is as follows: the experimental limits on $Br(h\to SS)$ are already strong enough to limit the annihilation rate of $SS$ pairs to the SM states to be much less than one picobarn, and consequently the metastable abundance of $S$ particles per nucleon is quite high, $Y_{p}\gg 1$. This leads to a massive injection of nucleons and mesons at early times, which raises the $n/p$ ratio, and creates larger yields of 4He compared to SBBN. Contributions of very light $S$ particles to the Hubble rate during the $n/p$ freeze out also raises $Y_{p}$. The limits on $\tau_{S}$ are robust, and have rather mild dependence on the uncertainties in our treatment. This is because the initial large metastable $Y_{S}$ abundance needs to be depleted prior to the $n/p$ freeze out time $t_{n/p}$, leading to the requirement $\tau_{S}\ll t_{n/p}$. Consequently $O(1)$ variations in the yields of mesons and nucleons in the final states can be compensated by small variations in $\tau_{S}$, parametrically on the order $\log^{-1}(t_{n/p}/\tau_{S})$, to produce the same influence on BBN. For the same reasons, our limits are also very insensitive to the exact observational constraint on $\Delta Y_{p}$, and we take a rather conservative limit of 0.01 (allowing $\pm 0.01$ deviations frrom the observed/calculated mean). From the point of the LHC physics, the most promising is a scenario with a mass $m_{S}$ not far below $m_{h}/2$. In that case, the effective decay length has to be on the order or smaller than $\sim 10^{8}$ meters, Fig. 6, providing a $10^{-6}$ minimum probability for a decay within a 100 m length purposely built detector. Given that the high-luminosity LHC would produce copious numbers of the Higgs bosons, there is a chance to cover the entire lifetime range for masses within $10-$to$-50$ GeV range.

It is easy to see that the above considerations can be generalized to other models of the Higgs-portal-coupled particles. For example, consider a fermion $\chi$, coupled to the Higgs via $H^{\dagger}H(\bar{\chi}\chi)$ or $H^{\dagger}H(\bar{\chi}i\gamma_{5}\chi)$ dimension-five operators, and having a small decay term such as e.g. neutrino portal $LH\chi$. The main analysis of our work can be recast for that model, especially in the part that connects Higgs decays with a metastable abundance of $\chi$. Evidently, for $Br(h\to\chi\bar{\chi})\sim Br(h\to SS)$ input, one will end up with $Y_{\chi}\sim Y_{S}$. The only change will be in the yields of mesons and baryons in the decays of $\chi$ compared to $S$. However, it is well known that already for $m_{\chi}$ above 250 MeV, the yields of pions and kaons is substantial Alekhin:2015byh , giving confidence that for the most parts same constraints we have derived for $\tau_{S}$ will translate to similar limits on $\tau_{\chi}$.

The analysis performed in this paper can be easily generalized to other models of metastable particles, with different types of interactions, via $Z$, $Z^{\prime}$ etc. In the limit when $Z^{\prime}$ is outside of the LHC reach, one could have a set of effective operators connecting $\chi$ with the SM fields, such as $\frac{1}{\Lambda^{2}}\bar{\chi}\gamma_{\mu}\chi\bar{q}\gamma_{\mu}q$, where $\Lambda$ is some energy scale. The $\chi$ pair-production cross section in this case will scale as $\sigma_{q\bar{q}\to\chi\bar{\chi}}\propto E_{q}^{2}\Lambda^{-4}$, where $E_{q}$ is a typical (anti-)quark energy, while cosmological annihilation cross section has $\sigma_{\chi\bar{\chi}\to q\bar{q}}v\propto m_{\chi}^{2}\Lambda^{-4}$ scaling. Therefore, the LHC-relevant cross section can be enhanced relative to the annihilation rate by a parametrically large ratio, $E_{q}^{2}/m_{\chi}^{2}$ if $m_{\chi}$ is parametrically smaller than the TeV scale. Therefore, one can easily have a range of parameters with a relatively large $\chi\bar{\chi}$ pair-production cross section, while having very small annihilation rates, rendering $Y_{\chi}\gg 1$, and resulting again in strong BBN constraints on lifetimes, $\tau_{\chi}<0.1$ seconds. Therefore, we conclude that some simple $Z^{\prime}$ mediated models of metastable particles can also be strongly restricted by cosmology, making them a perfect candidate for the searches of metastable particles at the LHC.

It is also instructive to consider models where constraints on the lifetime of metastable particles are much weaker. Clearly, one needs an effective new mechanism for the self-annihilation in the early Universe, as the Higgs channel is too inefficient. Staying within the Higgs portal models, consider the following potential with two real scalars,

[TABLE]

with the following hierarchy of couplings:

[TABLE]

These choices will lead to a long-lived $S_{1}$, somewhat shorter-lived $S_{2}$, a predominant decay of the Higgs boson to pairs of $S_{1}$, and cosmological depletion of $S_{1}$ via $S_{1}S_{1}\to S_{2}S_{2}$ annihilation with potentially a large cross section due to a sizeable $\lambda_{12}$ coupling. Most importantly, in this model the Higgs decay to pairs of $S_{1}$ does not result in a prediction of $Y_{S1}$ abundance, which can be quite small even for small values of $Br(H\to S_{1}S_{1})$. If $Y_{S1}\ll 1$, there would not be enough decay mesons and nucleons to affect early $n/p$ freeze out, and constraints on $\tau_{S1}$ will be coming only from the considerations of late decays with hadronic or electromagnetic energy injection. Instead of $\tau_{S}<0.1$ sec, one expects to have sensitivity to $\tau_{S1}\sim 10^{3}$ seconds, or even worse, beyond $10^{4}$ seconds, if decays of $S_{1}$ are mostly leptonic. This example is not unique, and there are other models where constraints on lifetimes and decay lengths are relatively lax, provided that there are extra channels that ensure efficient cosmological annihilation of metastable particles.

Acknowledgements

We thank D. Curtin, M. McCullough, P. Meade, M. Papucci and J. Shelton for soliciting this study, as well as D. Curtin and B. Shuve for very helpful discussions. Research at Perimeter Institute is supported by the Government of Canada through Industry Canada and by the Province of Ontario through the Ministry of Economic Development & Innovation.