Higgs phenomenology as a probe of sterile neutrinos

Jonathan M. Butterworth, ∗ Mikael Chala, 3, † Christoph Englert, ‡ Michael Spannowsky, § and Arsenii Titov ¶ Department of Physics & Astronomy, University College London, UK Institute for Particle Physics Phenomenology, Durham University, Durham DH1 3LE, UK CAFPE and Departamento de F́ısica Teórica y del Cosmos, Universidad de Granada, E–18071 Granada, Spain SUPA, School of Physics & Astronomy, University of Glasgow, Glasgow G12 8QQ, UK


I. INTRODUCTION
So far, no departure from Standard Model (SM) predictions has been established in collider experiments near or beyond the electroweak (EW) scale. This observation suggests that any new physics beyond the SM (BSM) is either very weakly interacting, or arises at a scale Λ much larger than the electroweak scale, or both. While the scenario with only new heavy physics is successfully described using an effective field theory (EFT) framework, the so-called SMEFT [1] (for a review see Ref. [2]), in the latter case, where new physics manifests itself in the presence of very heavy resonances outside the kinematic reach of the LHC on the one hand and light very weakly coupled degrees of freedom on the other, to describe the resulting BSM phenomenology the EFT framework involves not only SM fields but also new degrees of freedom, which are likely singlets under the SM gauge group. 1 One popular scenario is the EFT of the SM extended with a sterile neutrino N, also dubbed νSMEFT [5][6][7][8]. 2 Sterile neutrinos are present in many SM extensions, which aim to explain the origin of light neutrino masses. In particular, they are the main ingredient of the type-I seesaw [11][12][13][14][15] as well as the inverse [16][17][18] and linear [19][20][21] seesaw mechanisms. Although "canonical" (type-I) heavy neutrinos have masses close to the grand unification scale, mostly sterile neutrinos with much smaller masses can exist leading to various experimental signatures. Avariety of LHC studies have explored the phenomenology of the νSMEFT for N produced via contact interactions [5,[22][23][24], in W and top decays [24][25][26] as well as via Higgs decays with N decaying leptonically [27]; see also Ref. [28].
In this article, we focus on the production of one or two N via the Higgs with each N decaying into a photon and missing energy. We show that current data are not sensitive to the operators triggering these novel and clean Higgs signatures. Moreover, we go beyond the aforementioned works on this topic by performing much more realistic simulations of signals and background and therefore of the LHC reach. We also comment on the sensitivity that can in principle be gained using data-driven approaches in these clean final states.
This article is organized as follows. We introduce the νSMEFT in Sec. II and single out those operators which are not yet constrained by low-energy data. In Sec. III we discuss the potential of existing LHC searches and measurements to probe the aforementioned operators. Likewise, in Sec. IV we propose dedicated searches in monophoton and diphoton Higgs decays. We conclude in Sec. V.

II. FRAMEWORK
We consider the SM extended with one Majorana righthanded (RH) neutrino N and assume that its mass is below the scale of new physics Λ. The renormalizable Lagrangian gets modified as follows: where L SM stands for the SM Lagrangian, H represents the Higgs doublet, and L is the doublet of left-handed (LH) leptons with i ¼ e, μ, τ. Following standard notation, we have definedH ¼ iσ 2 H Ã and N c ¼ CN T with C being the charge-conjugation matrix. Also, m N is the Majorana mass of N.
Parametrizing new physics effects in terms of higherdimensional operators, at dimension five we have [6] where We note that for a single RH neutrino N the operator O NNB ¼ N c σ μν NB μν vanishes identically. At dimension six we consider the operators involving the Higgs doublet. The relevant Lagrangian reads where [8] and we have assumed the coefficients α to be real. In these equations, σ I with I ¼ 1, 2, 3 are the Pauli matrices, and For m N ≲ 10 keV, there are very stringent constraints on the new physics scale Λ from the cooling of red giant stars, implying Λ ≳ 4 × 10 6 TeV [6]. In the range 10 keV ≲ m N ≲ 10 MeV, supernovae cooling produced by the transitions νγ → N provides the strongest bound, which depends on m N as Λ ≳ 4 × 10 6 × ffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi m ν =m N p TeV [6], where m ν is the light neutrino mass. Taking m ν ∼ 0.01 eV, we find Λ ≳ 4 × 10 3 TeV (126 TeV) for m N ¼ 10 keV (10 MeV).
We are therefore interested in the regime in which N is relatively light but 0.01 GeV ≲ m N ≲ 10 GeV. The main decay channel is N → νγ induced by O NB and O NW [22]. The corresponding decay rate reads In principle, four-fermion operators are also present at dimension six, and can be expected to be more sizable because they can arise at tree level. However, there are models in which four-fermion operators are not generated at tree level; see the Appendix. Moreover, even if present these operators do not interfere with Higgs processes, which are the ones we are more interested in. Likewise, N → νγ is the dominant decay irrespective of the value of four-fermion interactions in the range of mass under consideration.
In this range of m N and assuming a standard cosmological history, the contribution of the new neutrino to N eff does not saturate the current Planck limit ΔN eff ≲ 0.3 [29]. (In alternative cosmologies, for example if the reheating temperature is close to ∼10 MeV, N does not achieve thermal equilibrium with the SM fields at any time, and cosmic microwave background constraints can be even more easily avoided [30]. ) We will make a number of assumptions on the coefficients of the considered operators to avoid constraints from low-energy data. First of all, we neglect the Yukawa couplings λ i in Eq. (1), which after electroweak symmetry breaking (EWSB) generate mixing of N with the SM neutrinos. Through the operators O i NB and O i NW this mixing would induce magnetic moments for the SM neutrinos, which are strongly constrained by reactor, accelerator and solar neutrino data [31,32]. This is entirely due to the missing t-channel mass suppression of the photon exchange. 3 Naively, even if λ i ¼ 0, the mixing would be induced after EWSB by the dimension-six operators O i LNH . However, without loss of generality, we can redefine the couplings λ i in Eq. (1) as λ i → λ i þ α i LNH v 2 =ð2Λ 2 Þ from the beginning, with hH 0 i ¼ v= ffiffi ffi 2 p being the Higgs vacuum expectation value (VEV). Such parametrization of the Yukawa couplings ensures that setting λ i ¼ 0 leads to no mixing, and we assume this in what follows. 4 Then, the Higgs-neutrino interaction reads Since the term of O i LNH proportional to v 3 cancels out, these operators contribute to neither the W AE → l AE N nor the Z → νN decay widths.
Upon EWSB, the operator O NNH contributes to the Majorana mass of N. Similarly to the discussion above, we can redefine m N in Eq. (1) In this way, m N is the physical mass of N. The hNN interaction arising from this operator has the form We attribute the smallness of ν i masses to the extremely small values of the coefficients α ij LH of the Weinberg operator given in Eq. (3). Therefore, we do not consider this operator in what follows.
In light of the previous discussion, we focus on the operators given in Eqs. (4) and (6)- (10). The operator O HN leads to the decay Z → NN with the following width (we neglect m N in analytical computations below): The operators O i NB and O i NW give rise to The operators O i HNe and O i NW contribute to the W decay width: It proves convenient for our further discussion to define the following operators: where c W ≡ cos θ W , s W ≡ sin θ W , with θ W being the Weinberg angle. Then we can rewrite Given these equations, we focus on a regime with α HN ¼ 0. This coefficient is extremely constrained by the measurement ΓðZ → ννγγÞ=Γ total (17)]. In this way we avoid the strong constraints on Z → γ þ p miss T [38][39][40][41]. Finally, we also assume that α i HNe ¼ 0, such that we completely avoid bounds from measurements of the W width. Moreover, the operators 3 On the other hand, massive mediators are largely unconstrained. Accelerator experiments (e.g., Ref. [33]) are relevant in a mass region of around less than 50 MeV, while masses in the keV range are subject to tight reactor data constraints (e.g., Ref. [34]). These are mass scales which are far below the typical hadron collider momentum transfers of Oð100 GeVÞ once trigger and selection criteria are included, which means that our results are insensitive to the concrete mass choice of neutrinos with masses ≲1 GeV. In this context, the low-energy measurements are only relevant when we make a concrete choice of small masses that do not impact our LHC analyses for the range that we consider. We will therefore not include the low-energy constraints explicitly in this work. 4 If the condition λ i ≈ 0 holds at a scale Λ ≫ v, then the different renormalization-group-equation running of the Yukawa and the O LNH operators might induce a SM neutrino dipole moment, O νA ¼ ðν iL σ μν ν c jL ÞA μν . The size of this operator can be estimated to be On the other hand, the latter is experimentally bounded to be α νA ≲ 2 × 10 −14 TeV −1 [31,35]. Therefore, λ i must vanish close to the EW scale (note that, below the EW scale, the dipole moment does not renormalize [36] O i HNe do not contribute to the processes we analyze in this work, so we can set α i HNe ¼ 0 without loss of generality. Under these assumptions we can express from Eq. (23) Thus, we are left with only α i NA , α i LNH and α NNH . The value of α i NA =Λ 2 is constrained by the measurement of Γ total W ¼ 2.085 AE 0.042 GeV to be ðα i NA =Λ 2 Þ ≲ 4π TeV −2 , so we can vary it in ½0.001; 4π for Λ ¼ 1 TeV. We have assumed that all three α i NA are of the same order. The lower bound is set by the requirement that N decays promptly enough (within 4 cm); see Eq. (13). Other low-energy constraints on α NA are very weak; see e.g., Ref. [42]. 5 For Λ ¼ 1 TeV, the coefficients α i LNH and α NNH can run in the ranges [0, 0.5] and [0, 0.05], respectively. The stringent upper bounds on these coefficients follow from the requirement that the partial decay widths of the Higgs boson in Eqs. (29) and (30) discussed later do not exceed the total Higgs width in the SM, Γ total H ≈ 4 MeV.

III. EXISTING SEARCHES
The most important processes at the LHC triggered by those operators which are very weakly constrained by low- (Note that the first and third processes do not interfere because the Higgs production is mostly initiated by gluons, while Drell-Yan production is initiated by quarks.) Let us first focus on the neutral-current Drell-Yan process.
In the limit of vanishing masses, the differential cross section for qq → ν i N is found to be where α ¼ e 2 =ð4πÞ is the fine-structure constant, and Q is the electric charge of the quark q. The integrated cross section is independent of s. For the process of interest we have where BðN → νγÞ ≈ 1 for the considered range of m N . This process can be constrained at the LHC in searches for events with one photon and missing energy. To the best of our knowledge, the most up-to-date search in this respect is the CMS analysis of Ref. [43]. Most importantly, this analysis requires exactly one photon with p γ T > 175 GeV and jη γ j < 1.44 as well as missing energy p miss T > 170 GeV. The ratio p γ T =p miss T is required to be smaller than 1.4 in order to reduce the background from γ þ jets. With the same aim, events are rejected if the minimum opening angle between p miss T and the transverse momentum of the four hardest jets is less than 0.5. (Only jets with p j T > 30 GeV and jη j j < 5 are considered in this cut.) Likewise, Δϕðp γ T ; p miss T Þ > 0.5. Finally, events are also rejected if they contain any electron or muon with p T > 10 GeV within ΔR > 0.5 from the photon.
The analysis considers two signal regions depending on whether j sin ϕj is smaller or larger than sin (0.5), which are further split into six p γ T bins in the range [175, 1000] GeV; see Table 1 in the experimental report [43].
We recast this search using dedicated routines based on ROOT v5 [44,45], HepMC v2 [46] and FastJet v3 [47]. Jets are built using the anti-k t algorithm with R ¼ 0.4, and defined by p j T > 10 GeV. For photons and leptons we require p T > 10 GeV. These objects are experimentally under very good control [48].
We find that the most constraining signal region is that with j sin ϕj < sinð0.5Þ and p γ T ∈ ½300; 400 GeV. The experimental collaboration reports the observation of 44 events, while 46.6 AE 4.0 are predicted in the SM. Using the CL s method [49,50], including the uncertainties in the estimation of the SM background, we obtain that the maximum number of signal events in this bin is 16. We estimate that the efficiency for selecting signal events in this bin in Drell-Yan processes triggered by O NA is ∼0.057. Using the leading-order (LO) production cross section before cuts, we obtain that α NA > 0.88 is excluded already at the 95% C.L., assuming that N couples to only one lepton family.
Interestingly, when running the simulated analysis over events of the type pp → h → ννγ, we find that none of the bins in this search constrain the operator O LNH . This can be easily understood because the distribution of the transverse which vanishes because h0jV μ jπ AE i ¼ f π p μ , with p being the pion four-momentum. momentum of the photon falls much faster in this process; see Fig. 1.
On another front, differential cross section measurements have been made at the LHC for a wide range of potentially relevant final states. These measurements are generally made in well-defined fiducial kinematic regions, giving them a high degree of model independence, and corrected for detector effects (such as resolution and reconstruction efficiency) within these regions, meaning they can be readily compared to generated signal events. We have used the CONTUR [51] tool to study whether our model would have had a visible impact on any of these measurements, all of which are currently consistent with the SM. To do this, we use the Herwig 7.1.5 [52,53] event generator to read the UFO [54] files of our model and produce simulated collision events. All processes with any BSM content in the final state, or on-shell intermediate states, are generated at LO tree level. For each parameter point, one million events are generated, implying an integrated luminosity at least equivalent to that of the data, and typically much higher. These events are then passed to RIVET [55], which contains implementations of a large number of the relevant LHC analyses as well as the measurement data derived from HEPDATA [56]. The effect of injecting signal events on top of the data for all these measurements is then evaluated by CONTUR, and the most sensitive distributions for any given model point are identified and used to derive a potential exclusion, using a χ 2 test.
Scanning the range 10 −3 < α NA < 4π with the other couplings set to zero, we observe that for α NA ≳ 3, heavy neutrino production via Drell-Yan is significant, q 0q → Ne AE and qq → Nν e . For example, at α NA ¼ 3.4 and m N ¼ 0.2 GeV, the inclusive cross section for these processes combined is 40 pb, dominated by the channels involving e AE . As they contain an electron and missing transverse energy, these events can easily populate the fiducial phase space of measurements aimed at W bosons decaying to electrons and neutrinos [57], with the photon from the N decay also meaning that they can impact upon W þ γ measurements [58]. They, and the neutral-current Drell-Yan processes, can also contribute to inclusive photon and photon-plus-jet measurements [59,60], where no veto is made on the rest of the event, and photon-plus-missingenergy measurements [61]. The resultant exclusion is shown in Fig. 2. More recent measurements, and those from CMS in these channels, are not yet included in RIVET and so are not used.
We note that there is no dependence on m N over the range considered. Setting α NA ¼ 1, α NNH ¼ 0 and scanning 0 < α LNH < 0.5 for the same range of m N , some events do enter the fiducial region of the same measurements, but there is no impact at the 1σ level or above. The same is true for scanning 0 < α NNH < 0.05 with α NA ¼ 1, α LNH ¼ 0. The fact that some events do populate the acceptance of these measurements means there may be sensitivity as the measurement precision is increased with higher integrated luminosity.
These existing limits leave a large part of well-motivated parameter space uncovered. To probe this region, new analyses, not yet conceived, are required. 6 We discuss such analyses in the next section.

IV. HIGGS SEARCHES IN THE MONO AND DIPHOTON + MISSING ENERGY CHANNELS
In this section we investigate the sensitivity reach of the LHC to the mono-and diphoton þ missing energy channels through novel Higgs decays.
There are different strategies to constrain new physics signals at colliders. On the one hand, if a good understanding of the background and the signal can be achieved this can be used to inform an experimental search in cutand-count or more sophisticated multivariate analyses, in line with the previous section. This approach is the major driving force behind searches in complicated multiscale final states. As also discussed in the previous section, precision differential measurements of relatively simple final states can also contribute.
On the other hand, if a signal process is clustered in a particular phase space region (e.g., for resonance searches) we can use sidebands to constrain the background with a minimum of theoretical input using data-driven approaches. Such a strategy comes into its own when the final-state objects are experimentally under good statistical and systematic control, which is the case for photons already at low transverse momenta [48]. If the shape of a background can be accurately described using fitting techniques across a large range of a kinematical observable, We have generated events setting α NA ¼ 1 for both processes, and in addition α LNH ¼ 0.1 for the second process. All the other coefficients have been set to zero. 6 The very recent search for Zh with h → γ þ p miss T by the CMS Collaboration [62] had very little sensitivity. such a strategy can be used to detect small signals on top of large backgrounds even when the latter are theoretically not well understood. A prime example of this strategy is arguably the Higgs discovery via its decay to γγ with relatively small signal-to-noise ratio [63,64].
Our search, although involving missing energy with different systematic properties, shares many similarities with the h → γγ search: the background is large and the expected branching ratio h → γðþγÞ þ p miss T is small, i.e., the new physics signal is likely to be dwarfed by the expected theoretical uncertainties associated with monoand diphoton production in a complicated hadronic environment. However, the signal normalization is accurately known (see e.g., Refs. [65][66][67][68]) and its relevant final-state kinematics is entirely determined by the Higgs mass which can be accurately extracted from subsidiary measurements; see for example the ATLAS and CMS combination of Ref. [69]). This implies a distinct shape of the new physics signal, i.e., there is a Jacobian peak in the transverse momentum distribution of the photon or photon pair, depending on which final state we are interested in as shown in Fig. 3. Note that this way the resonance cross section extraction is also not impacted by BSM contributions to the continuum as indicated in Fig. 1: the background shape might change but the resonance will still have a distinct shape, which can be extracted from the continuum for large enough data sets.
In the following we take inspiration from the h → γγ search and estimate the sensitivity at the high-luminosity LHC by performing a template fit on the signal and background distributions using RooStats [70] (with background shape estimates taken from Monte Carlo), leaving their normalizations as free parameters. (Using the same approach, we are able to reproduce e.g., the expected FIG. 2. CONTUR exclusion in the α NA , m N plane, for Λ ¼ 1 TeV. The left-hand inset shows the 2-and 1-σ exclusion contours based on the heat map on the right. Within the diagonal cyan region the heavy neutrino would decay within the detector volume; above, it is effectively prompt and below, it is effectively stable. ATLAS p-value of the 8 TeV h → γγ search of Ref. [63] within 10%. This highlights that such an approach is highly feasible when all experimental aspects are under good systematic control, which we assume here implicitly, but not unrealistically.) The 95% C.L. constraint on the signal modifier when agreement with the background-only hypothesis is given can then be understood as a direct constraint on the respective branching ratios when using the signal normalization of pp → h from Ref. [68]. The expected background cross sections for the inclusive selection criteria that underpin Fig. 3 are σðγ þ p miss T Þ ≃ 14 pb and σðγγ þ p miss T Þ ≃ 10 fb. This way we obtain at 3 ab −1 using signal and background templates generated with MADGRAPH [71] as shown in Fig. 3. In line with the previous section we require a minimum p T of the photon of 10 GeV for our mock γ þ p miss T data and the expected SM h → Zγ contribution is subtracted from these numbers. Note that owing to the much smaller expected background of the γγ þ p miss T analysis, the signal is naively easier to isolate, but it is considerably more washed out due to missing energy systematics. These issues fall into the area of experimental expertise, and hence we limit ourselves to the sensitivity estimated along the lines above, but we choose harder photons p T ≥ 15 GeV with separation in the pseudorapidity-azimuthal angle plane of at least 0.4 as well as a minimum missing transverse energy of 25 GeV for the diphoton analysis to partially take into account the more complicated nature of this process and report results separately.
The bounds in Eqs. (27) and (28) can be translated to constraints on α LNH , α NNH and α NA by means of (See Ref. [72] for former partial computations.) Each of the three Wilson coefficients can be bounded independently by setting the remaining two to zero. (Note that any other choice would lead to a more stringent constraint.) The sensitivity at the high-luminosity LHC is listed in Table I. We remind the reader that all these prospects apply only if α NA =Λ 2 ≳ 0.001-0.1 TeV 2 , depending on m N ; see Fig. 2.

V. CONCLUSIONS
In summary, we have studied the phenomenology of the SMEFT extended with a light RH neutrino N in the regime in which the latter decays almost exclusively into a photon and a neutrino.
Using low-energy and LHC data such as measurements of the W, Z and Higgs bosons, bounds on neutrino dipole moments, measurements of SM differential distributions at the LHC (as implemented in CONTUR [51]), as well as searches for single photons with missing energy [43], we have singled out those directions not yet constrained. They include mostly operators triggering new Higgs decays, namely h → γ þ p miss T and h → γγ þ p miss T . We have subsequently provided new search strategies to be performed at the LHC sensitive to the aforementioned processes. For order-one couplings, we have shown that, with 3 ab −1 of data, these analyses can potentially unravel new physics at scales Λ ≲ 2 TeV (2000 TeV) for leptonnumber-conserving (-violating) operators. For comparison, let us add that searches for h → NN triggered by O NNH , with N → qql are expected to test scales as large as ∼100 TeV [27]. Likewise, top decays into blN, mediated by four-fermion operators and with long-lived N have been shown to probe only Λ ≲ 1 TeV [24].

APPENDIX: MODEL
Let us consider the SM extended with two vectorlike fermions X E ∼ ð1; 2Þ 1=2 , X N ∼ ð1; 1Þ 1 and a singly charged scalar φ ∼ ð1; 1Þ 1 . The numbers in parentheses and the subindex represent the quantum numbers under ðSUð3Þ c ; SUð2Þ L Þ and the hypercharge, respectively. We also assume that these new fields are odd under a Z 2 symmetry under which all SM fields as well as N are even.
The new relevant Lagrangian reads Let us focus on the regime M X E ; M X N ; M φ ∼ M ≫ v, g N ≪ g L ; g X . The new particles can be integrated out before EWSB by matching (off-shell) amplitudes in the UV to the corresponding amplitudes in the EFT. One can easily check that tree-level operators vanish.
Therefore, we concentrate first on the amplitude given by the diagrams in Fig. 4. Using p 2 , p 3 and p 4 as independent four-momenta (p 1 ¼ p 2 þ p 3 þ p 4 ), and to first order in p i we get and A n , B n and C n are four-dimensional integrals defined by (see e.g., Ref. [73])  k μ k ν k ρ k σ ðk 2 − M 2 Þ n ¼ 1 4 C n ðg μν g ρσ þ g μρ g νσ þ g μσ g νρ Þ: ðA6Þ Explicitly, Γðn − 3Þ ΓðnÞ ; ðA8Þ Adding all pieces together and simplifying further, we get The only operator in the IR that contributes directly to the same amplitude is O NA ; it reads Upon requiring M UV ¼ M IR we finally obtain α NA Λ 2 ¼ g L g X g N e 192π 2 M 2 : In order to obtain α LNH , we compute the amplitude given by the diagrams in Fig. 5 to zero momentum. We have iM 0 a ¼ − 3 ffiffi ffi 2 p g N g 3 X g Lū ðp 5 ÞP L ½6C 5 þ 24M 2 B 5 þ M 4 A 5 uðp 1 Þ; ðA13Þ ffiffi ffi 2 p λ φH g X g L g Nū ðp 5 ÞP L ½4B 4 þ M 2 A 4 uðp 1 Þ: Adding both contributions, we obtain ffiffi ffi 2 p π 2 M 2 ½λ φH − g 2 X ūðp 5 ÞP L uðp 1 Þ: ðA15Þ In the EFT we obtain instead from which we obtain Redundant operators can be generated in the off-shell matching, therefore potentially contributing to O NA and O LNH after using the equations of motion of the SM þ N. The relevant list of such operators reads  ðA21Þ (The addition of Hermitian conjugates is implied.) Other operators (not related to the ones above by algebraic identities or integration by parts) involve two copies of N, as we have cross-checked using BasisGen [74]. Therefore, their Wilson coefficients are suppressed by two powers of g N and therefore negligible within our approximation g N ≪ g L ; g X . We expect contributions to O NA to be small. Likewise, further contributions to α LNH come from the equations of motion of the Higgs [1]: They are therefore suppressed by a further factor of λ H ∼ 0.1 and therefore negligible. In summary, assuming Λ ¼ M, in good approximation α NA ∼ g L g X g N e 192π 2 ; α LNH ∼ g N g L g X 96π 2 ½λ φH − g 2 X : ðA23Þ We note that in the strongly coupled regime, and for m N ∼ 1 GeV and Λ ¼ 1 TeV, N decays effectively promptly within the detector (see Fig. 2), and α LNH is within the reach of our analyses (see Table I). For example, neglecting λ φH and for g L ¼ g X ¼ ffiffiffiffiffi ffi 4π p and g N ¼ 1, we get α NA ∼ 0.002 and α LNH ∼ −0.17. α NA grows up to ∼0.03 for g L ¼ g X ¼ 4π.