Probing ALP-Sterile Neutrino Couplings at the LHC

In this work, prospects to probe an overlooked facet of axion-like particles (ALPs) -- their couplings to sterile neutrinos -- are presented. We found that mono-photon searches have the potential to constrain ALP couplings to sterile neutrinos when a new heavy scalar boosts the ALP decay yields. Working within an effective field theory (EFT) approach, we scan the parameters space to establish the reach of the 13 TeV LHC to probe such couplings. We found that regions of the parameters space evading several experimental constraints that can be probed at the LHC. Moreover, a complementary role between the LHC and various experiments that search for axions and ALPs can be anticipated for models where ALPs interact with sterile neutrinos. We also present the UV realization of a model having an axion-like particle, a heavy scalar and sterile neutrinos whose parameters are spanned by our EFT approach. The proposed model contains a type of seesaw mechanism for generating masses for the active neutrinos along with sterile neutrinos involving the high energy scale of the spontaneous breaking of the global symmetry associated to the ALP. Some benchmark points of this model can be discovered at the 13 TeV LHC with 300 fb$^{-1}$.


INTRODUCTION
The study of pseudo-Nambu-Goldstone bosons such as axions and axion-like particles (ALPs) has lead to many developments in high energy physics. Theoretically, these particles are expected to have two main features: they are light in comparison to the scale of spontaneous breaking of the global symmetry which originates them; and their couplings to the Standard Model (SM) particles are, at least, suppressed by the inverse of that same scale. A well known solution to the strong CP problem is the Peccei-Quinn mechanism which is based on the implementation of an anomalous chiral global U(1) P Q symmetry [1,2], whose spontaneous breaking gives rise to the axion [3,4].
The axion gets mass inversely proportional to the scale of spontaneous breaking of the U(1) P Q symmetry due the peculiar feature that this symmetry is also broken explicitly at the quantum level.
The ALP has some interactions similar to the axion but, differently, its mass is not directly related to its couplings with the SM particles. In the ALP case, the associate spontaneously broken global U(1) symmetry may be generic, and also approximate due the existence of operators breaking it explicitly. Such operators give to the ALP a mass that is not necessarily related to the scale of the spontaneous broken of the U(1) symmetry. Even in this situation, the ALP mass can be considered naturally small by the argument that in the limit where they vanish, the symmetries are augmented [5]. The search for axions and ALPs is currently among the main experimental efforts in particle physics. For the general aspects of axions and ALPs see [6][7][8][9][10], for example. In this work we will deal solely with ALPs.
Many experiments and astrophysical observations furnish constraints on ALP masses and couplings to photons, gluons and electrons [6]. However, the couplings of ALPs to neutrinos and other dark particles have been overlooked in comparison to their couplings to more easily observable probes like photons and electrons. This is comprehensible though. In colliders, the search for a light scalar decaying invisibly is difficult in the gluon fusion production mode once the missing energy spectrum produced along with a visible jet, in a monojet-type signature, is expected to be soft. A good probe for the ALP(a)-photon coupling is Z → γ + a → 2(3)γ [11,12] but if a decays invisibly, the monophoton signature would produce a not too hard photon that would be difficult to observe in colliders. In this work, we not only explore a mass region that is more promising for collider searches [13] but also consider a class of models where particles present harder spectra.
Several studies have already considered the possible signals and prospects of ALPs in colliders for specific couplings with gauge bosons and fermions. For example, depending on the scenario considered, ALPs can be constrained through: Z boson decay into two or three photons [11,12] and other SM heavy resonances decay [14,15]; monophoton, monojets, and triphotons [16], decay into two photons and two leptons [13], non-resonant production [17,18], and photon pair production in ultra-peripheral collisions [19] signals at the LHC; production via Primakoff process in fixed target experiments [20][21][22] such as NA62 and SHiP [23,24], in the FASER experiment [25] and also in photon-beam experiments [26]; flavor violation interactions via couplings to the W ± bosons [27,28] and leptons [29]; ALPs as mediators of interactions between the SM particles and dark matter [30].
Realistic models embodying ALPs predict at least one-loop couplings to all the Standard Model particles. At the tree level, ALPs are expected to participate in the scalar potential. Possible ultra-violet complete theories involving ALPs, by their turn, include new heavy fermions which might need new scalars if their masses are generated via a sort of spontaneous symmetry breaking mechanism. In those scenarios, a large scalar sector with the pseudoscalar ALPs and another new heavy scalar, at least, is plausible. Augmented fermionic and scalar sectors are very compelling sources of new physics which might be necessary in order to stabilize the SM potential and/or to trigger baryogenesis, for example.
Hidden valley, MSSM and NMSSM [31][32][33][34] are examples of models where a heavy scalar is predicted to decay resonantly into lighter scalars. For example, in NMSSM, a heavy CP -even scalar can decay into two CP -odd ones giving very distinctive signatures [35]. In the context of axion-like particles, a heavy scalar might decay into two ALPs which, if not long-lived, decay inside the LHC detectors.
The MiniBoone [36] and LSND [37] experiments have reported an excess of appearance of electron neutrinos in a muon neutrino beam for over a decade which, in combination, amount to a 6σ statistical significance level. Explaining this observation in a two-neutrino oscillation model would require, at least, four neutrinos, thus evidencing the existence of sterile neutrinos.
This might considered, so far, one of the most robust experimental evidence for physics beyond the SM paradigm. Recently, models where a sterile neutrino decays to a light scalar and SM neutrinos were proposed [38,39] in order to alleviate the tension between the MiniBoone/LSND explanation of those observed excesses and other experiments that have not been observed an excess of disappearance of muon neutrinos. This light boson could be an ALP.
In view of these motivations, we want to explore the possibility of a new channel not previously studied, the resonant ALP pair production where one of the ALPs decays to sterile neutrinos and the other one to photons. By sterile neutrino, we mean a generic long-lived neutral fermion which originates from a singlet field under the Standard Model gauge group.
We focus on models where a light ALP is allowed to decay to photons, sterile neutrinos, light charged leptons and hadrons. ALPs of order of GeV and heavier can be probed in gluon and photon fusion processes [40]. Decays to electrons open up interesting new signatures but we will not explore them right now. If the scalar that decays into ALPs is heavy and the ALP is short lived, they get boosted in the laboratory frame and decay inside the LHC detectors leading to a pair of collimated photons mimicking a monophoton signal -a photon jet. This channel has been intensely studied by the LHC collaborations in the search for dark particles [41], and also as signals of a high-mass resonance decaying into pairs of such effective photon jets via light scalar resonances [42].
Our goal in this work is to explore the 13 TeV LHC reach to probe ALP-neutrino couplings in the monophoton channel where one heavy scalar, S, decays to two light ALPs, aa, which then decay to photons and sterile neutrinos that pass the detector leaving with no trace. From now on, we assume that the ALPs are boosted and lead to very colimated photons which hit a single detector cell mimicking a monophoton channel. This kind of signature with colimated photons has been explored in previous works like [11,12].
We will work within an effective field theory formulation of the interactions between the ALPs, the new scalar S, the sterile neutrino N and the SM fields. However, we also provide an ultraviolet complete construction to exemplify a concrete model which motivates the choice of the EFT parameters. In this model, the heavy scalar and the ALP are the components of a same SM singlet complex field with self-interactions dictated by a scalar potential. It is reasonable to expect that the high energy scale in which the approximate U(1) symmetry related to the ALP is spontaneous broken participates in the neutrinos mass generation mechanism. Thus, the model is also built to furnish a mass generation mechanism that can lead to sterile neutrinos with masses at the MeV scale as well as active neutrinos with masses at the sub-eV scale. Moreover, we are going to assume that the same new heavy states mediate both the S and ALP 1-loop interactions to SM particles.
After integrating out these heavy states, the new physics scale Λ parametrizes all the new scalar interactions making our EFT model more predictive.
Our paper is organized as follows. In Section 2 we present the effective Lagrangian on which our developments are based and discuss the main parameters entering in the decay widths. In the Section 3 we discuss the implementation of the several existing experimental constraints. Our analysis on the search for the ALP-sterile neutrino interactions with monophotons at the LHC are presented in Section 4. A generic model for a complete theory at the ultra-violet including a seesaw mechanism to generate sterile neutrinos masses at the MeV scale as well as active neutrinos masses at the sub-eV scale is present in section 5. Our conclusions are given in Section 6.

EFT LAGRANGIAN AND PARTIAL DECAY WIDTHS
We parametrize the relevant interactions of the heavy scalar, S, the ALP, a, the sterile neutrino, N , and SM gauge bosons, electrons and muons, with the following effective Lagrangian We assume that all the scalar couplings are mediated by the same new physics scale Λ as discussed in the Introduction. The parameters c BB , c W W , c GG , k BB , k W W , k GG , c aN , c ae and c aµ are all dimensionless whereas f S has dimension of mass. The coupling to gauge bosons preserve the SM group symmetry. The field strength tensors of the U (1) Y , SU (2) L and SU (3) C gauge bosons are, respectively, B µν , W µν and G a µν . The dual field strength of the U (1) Y gauge boson field is defined asB µν = 1 2 αβµν B αβ , with the totally anti-symmetric tensor being such that 0123 = 1, and analogously to the other fields. In the Section 5 we give an example of a concrete ultra-violet complete model containing heavy fields which can lead to the low energy effective Lagrangian like the Eq. (2).
Using the equation of motion for the neutrinos, we see that the interaction term with them is actually proportional to their masses. That is the reason why we do not expect that couplings to the SM neutrinos can be probed at all, only to a heavy sterile neutrino of mass m N . An invisible ALP decay is therefore a clear signal of interactions with sterile neutrinos or some new dark state.
By its turn, the Saa coupling, f S , has dimension of energy and can, in principle, be large.
In spite of the simplifying assumption of a same new physics scale mediating both S and ALP interactions, the couplings of these new states to these scalars are not necessarily the same. Both S and a couplings to the SM particles are supposedly given in terms of the structure and parameters of a complete high energy theory. For example, if S and a are originated from scalar singlets under the SM group they might have tree level couplings with the known fermions and the Higgs boson proportional to mixing angles. In this case, their couplings to the vector bosons, γ, Z 0 , and W ± can arise through one-loop order in perturbation theory involving new particles carrying SM quantum numbers. In special, the ALP couplings depend in general on the associated U(1) symmetry charges of the fermions and their properties under the SM group (see Section 5). If S and a are the real and the imaginary components of the same complex scalar field, as the model presented in Section 5, the parameters c BB , c W W , c GG and k BB , k W W , k GG above can be related.
Anyway, we stick to the general case where these parameters can be all different from each other.
The partial widths of S and a are given right below In the formulae of Eqs.(3) above, m S and m a denote the heavy scalar and the ALP masses, respectively; m Z and m W the masses of the Z and W bosons, respectively; = e, µ; and s 2 w ≡ sin 2 θ W ≈ 0.231 is the sine squared of the electroweak mixing angle. In the ALP mass range from 0.4 GeV up to 2 GeV, the effective ALP-gluon coupling cannot describe the coupling to hadrons correctly as it does not take the mass thresholds into account. For this mass range we take the main hadronic modes individually in our calculations as showed in Ref. [43]. Yet, as we are going to see, after the neutral pion channels is open, the monophoton signature vanishes and the prospects for the LHC observation die out if m a 0.4 GeV effectively. Nevertheless, it is possible to choose a fermion content in a UV-complete realization of the model in which the couplings to gluons are suppressed raising the branching ratio of the other channels (see Section 5).
Three main factors dictate the number of monophoton signal events that can be produced in pp collisions at the LHC. The resonant production cross section of the ALPs, σ(pp → S → aa), the branching ratio BR(S → aa), and the product BR(a → γγ) × BR(a → NN ).
In Fig. (1), we display the branching ratios of the ALP for masses up to 1 GeV in eight different scenarios. In all of them, we fix m S = 1000 GeV, f S = 100 GeV, and k BB = k W W = k GG = 1.
The sterile neutrino mass is fixed to 0.1, 1, 10 and 100 MeV at the first, second, third and fourth row, respectively, and ALP-leptons coupling is set to 0 in the plots of the first column, and at 10 Showing scenarios for different ALP-leptons coupling anticipates the two frameworks that we are working with -the EFT approach, where c a is free, and the UV complete model of Section 5 where this coupling is suppressed. We basically see from Fig. (1) that switching off the leptons couplings leaves more room for γγ and N N decays naturally. In special, we are interested in the product BR(a → γγ) × BR(a → N N ), depicted as dashed black lines. If the sterile neutrino mass is larger than 0.1 MeV, we see that including leptons decays have a minor impact actually. This is because the ALP-sterile neutrino mass is proportional to its mass.
As m N , the sterile neutrino mass, increases, the region of ALP masses where decays to γγ +N N shrinks and as soon as m N approaches the hadrons threshold, the hadrons channel dominates and decays to photons and neutrinos become rare. Finally, notice that for m N < 100 MeV, there exist ALP mass regions where BR(a → γγ) × BR(a → N N ) can reach its maximum of ∼ 0.25. These spots will be the most promising ones to search for monophoton events at the LHC.
The dominant branching ratios of the scalar S, as a function of its mass, are shown in Figure (  Next, let us start discussing the search for signals of the model at the LHC by first evaluating possible constraints to the model.

EXPERIMENTAL CONSTRAINTS
The interactions of Eq. (2) can be probed in a number of ways and several experimental constraints apply: (1) Resonant production of the scalar S and decay to dijets, diphotons, ZZ, W W and Zγ. The monophoton search is effective as long as the ALP decay to hadrons is closed. In this case, dijets and monojet signals do not receive contributions from S → aa but, diphotons and monophotons do. Note that imposing constraints on the ALP-gluons coupling would be beneficial to monophoton searches.
(3) ALPs should be sufficiently short lived in order to decay inside the detectors.
(4) Constraints from Z decays to photons.
We discuss the details of each one of these constraints right below. All the following experimental constraints were obtained from the 13 TeV LHC data, except for the last one of Z boson decays to photons that come from the 8 TeV LHC run and Tevatron data.
1. Diphoton searches -The diphoton final state provides a clean experimental signature with excellent invariant-mass resolution and moderate backgrounds. According to [44], searches for new phenomena in high-mass diphoton final states with the ATLAS experiment at the LHC gives an upper limit on at 95 % confidence level.

Dijet searches -The search for resonances decaying into pairs of jets by the CMS [45], with
an integrated luminosity of 36 fb −1 , places the limit on the resonance production cross section times the branching ratio Br(S → jj) times the acceptance A given by at 95% confidence level. Here, acc is the acceptance for narrow resonances, with the kinematic requirements |∆η jj | < 1.3 for the dijet system and |η j | < 2.5 for each of the identified jets, where η j is the rapidity of a jet.
3. ZZ pairs -A search for a new scalar resonance decaying to a pair of Z bosons by the CMS [46] with an integrated luminosity of 35.9 f b −1 , where the Z boson pair decays are reconstructed using the 4l, 2l2q, and 2l2ν final states, imposes the following upper limit on the resonance production cross section times the branching ratio Br(S → ZZ) at 95% confidence level.
4. W W pairs -Searches for neutral heavy resonances performed in the W W → eνµν decay channel either directly or via leptonic tau decays with additional neutrinos, with an integrated luminosity of 36.1 f b −1 , were made by the ATLAS experiment at the LHC in [47].
According to these searches we must have the upper limit at 95% confidence level.
The search in the eνµν decay mode is complementary to searches performed in other decay modes. In particular, the sensitivity to low mass resonances is higher in the fully leptonic final state than in final states that include jets due to background from jet production.

5.
Zγ searches -The Z(→ + − )γ final state can be reconstructed with high efficiency and good invariant mass resolution. Based on Zγ searches of a heavy resonance by ATLAS [48], with 36.1 fb −1 , we set an upper limit on the resonance production cross section times the branching ratio Br(S → Zγ) given by at 95% confidence level.
6. Monophoton searches -Models that predict a photon with high momentum and large E T in pp collisions have been constrained by the LHC. The results of a search of such events in the ATLAS study of Ref. [49] with an integrated luminosity of 36.1 fb −1 gives us an upper limit on the total cross section for the process at 95% confidence level. Here, acc is the acceptance for monophoton final states with the kinematic requirements E γ T > 150 GeV, |η| < 1.37 or 1.52 < |η| < 2.37 and ∆φ(γ, E T ) > 0.4. Also it is required events without jets or just one jet with p T > 30GeV , |η| < 4.5 and ∆φ(jet, E T ) > 0.4, no charged leptons in the final state and E T > 300 GeV.
7. ALP lifetime -We require that the ALP decays inside the electromagnetic calorimeter which, at ATLAS and CMS, is situated approximately at 1 meter from the interaction point. The distance that the boosted ALP produced in the decay of the heavy scalar S travels from the interaction point before decaying is given approximately by [11,12] The total width of the ALP can be computed from where h i denotes a hadronic final state of total mass m h i . Substituting this equation in Eq. (10) leads to the constraint 8. Photonic decays of the Z boson -Important constraints from searches for Z boson decays into photons are given by the CDF Collaboration [50]: Br(Z → γγ) < 1.45 × 10 −5 and the ATLAS Collaboration [51]: The azimuthal opening angle between the two photons coming from the ALP in the decay Taking the resolution of the CDF detector to be 120 mrad [52], we see that the two photons from the ALP will be colimated if m a < 2.7 GeV, and the CDF constraints apply [12]. The resolution of the LHC detectors is ∆φ = 20 mrad. The search for 3γ decays of the Z boson performed by the ATLAS Collaboration [51] thus apply if the photons from the ALP decay are not collimated, that is it, for m a > 0.46 GeV.
In both cases, the relevant branching ratio of this Z decay channel is Neglecting the ALP mass in comparison to the Z boson mass we now have another constraint on the parameters in the EFT Lagrangian where Br exp Z is the CDF limit if m a < 2.7 GeV, and the ATLAS limit if m a > 0.46 GeV. As we are going to see, the CDF limit is, effectively, the bound to be respected. Now, we are going to investigate the prospects of the LHC to observe monophoton events from the production of ALPs taking into account all the constraints just described.

Sampling of parameters and ALP lifetime
We are interested in the regime where one ALP from S → aa decays to two colimated photons and the other one to neutrinos, a monophoton signal therefore, as we discussed earlier.
The EFT Lagrangian of Eq. (2) encompasses 11 independent parameters plus the S, N and ALP masses, totaling 14 parameters. Instead of choosing a few representative benchmark points, we performed a thorough exploration of the parameters space looking for insights of the interplay among those parameters to find the most promising spots for a collider search. In order to take into account all the constraints simultaneously, we generated 2 × 10 6 random sets of parameters and tested them against all the constraints described in Section (3). The range of the parameters drawn from uniform or log-uniform distributions are the following The range of the c and k coefficients are chosen with the assumption that they should lie in the perturbative regime of a concrete model. Also, we require c aN Λ m N < 1 in order to keep ourselves in the perturbative regime as well.
The ALP mass and its coupling to photons are very constrained by a variety of experiments.
In Ref. [13], a mass range of a few MeV is established as the one which is more promising for collider searches. The ALP-neutrino and ALP-leptons couplings, c aN , c ae , and c aµ , respectively, The requirement that the ALP must decay inside the detector, given by Eqs. (10) and (12), strongly constrains the ALP mass. Actually, only 30% of all points generated pass this constraint.
In Fig. (4), we see that the bulk of viable points that pass all the constraints lie above m a ∼ 100 MeV and m N ∼ 1 MeV. Notice that the neutrinos decay channel make it easier to respect this requirement once the bulk of the points lie in the region where m a > 2m N . Concerning the neutrino masses, considering only one type of neutrino, the branching ratio of the ALP into neutrinos turns out to be proportional to (m N /m a ) 2 . If m N << m a , Br(a → N N ) gets diminished, making the cross section too small to give any detectable signal. The cosmological upper limit on the sum of the SM neutrinos masses [6], ν m ν < 0.72 eV, by its turn, imposes an upper limit on the SM neutrino masses which make their branching ratios very small. This is an immediate consequence of the ALP-fermion interaction that scales with the fermion mass. For this reason, a missing energy signal at colliders can only come from ALP decays to a new heavier neutrino or dark fermion.
There are laboratory and cosmological constraints over the sterile neutrinos mixing with the active ones for the mass scales below GeV which we are considering here, as can be seen in Ref. [53] and references therein. We discuss these constraints in Section 5 where an ultraviolet completed model that can account for the effective interactions of the ALP and mass generation for the neutrinos is presented.

Details of the signal and backgrounds simulations
In order to estimate the reach of the LHC to constrain signals of ALPs in the monphoton channel, we simulated events for the signal pp → S → aa → γγ + N N , and its main SM backgrounds: (1) the irreducible Z(→ νν) + γ, and (2) the reducible W (→ lν) + γ with a missing electron or muon coming from the W boson, and (3) the reducible W (→ eν e ) where the electron is misidentified as a photon, following the ATLAS study of Ref. [54]. The simulations were performed at Leading Order using MadGraph5 [55]. We simulated parton hadronization and showering with In the subsequent analysis, we compute the number of signal events N S for an integrated luminosity L as follows  particles impact the cut efficiency and, as we discussed above, is affected only by m S .
The simulated events were generated requiring to pass the following basic cuts where p γ T and η γ denote the transverse momentum and pseudorapidity of the hardest photon of the event, respectively. In the case of the W → e + ν e background, we required p e T > 165 GeV and |η e | < 2.5. In the subsequent analysis, the electron from this background is treated as a photon as discussed next. Besides the main contributions of Zγ and W γ, the backgrounds might receive various subleading contributions. The most important comes from W → e + ν e where the hard electron is misidentified as a photon in the detector. The experimental study of Ref. [54] depicts all the subleading contributions, many of them are hard to be reproduced like, for example, spikes in the electromagnetic calorimeter or effects from the beam halo. We choose to rescale our number of backgrounds by the factor corr = 1 + subleading leading , where "subleading" is the total number of subleading backgrounds, and "leading" the number of Zγ and W γ events, respectively. The subleading contributions amount to around 30% of the total backgrounds where the electron misidentifcation contribution from W decays account for around half of the subleading contributions for the basic cuts of Eq. (18).
The additional cuts shown right below, still following Ref. [54], are imposed to clean up background events further where ∆φ( p T , p T ) is the angle between the missing momentum vector p T and the momentum vector p γ T of the photon and p j T of the leading jet in the transverse plane. Apart from the missing E T , all the other cuts are identical to those from Ref. [54].
The cut efficiencies and the cross section of production of the signal and backgrounds before and after cuts of Eq. (19) are shown in the Table 1 for E T > 170 GeV and E T > 500 GeV. The backgrounds are efficiently eliminated by hard missing transverse momentum cuts while the signals maintain a good cut efficiency especially if the S scalar is heavy.
The backgrounds normalizations after the cuts of Eq. (19) were taken from the ATLAS simulation of Ref. [54]. Our simulations for the Z + γ and W + γ backgrounds were used to calculate the correction factors for the cut efficiency in terms of harder E T thresholds keeping all the other cuts of Eq. (19). We checked that our normalizations after cuts, with the NLO QCD K-factors of 1.46 [55] and 2.41 [55] for the Z + γ and W + γ, respectively, are 0.7 and 2.7 times the ATLAS normalizations, respectively. The W → e + ν e background is severely affected by hard missing transverse energy cuts, nevertheless we assume a conservative approach concerning the subleading backgrounds and just rescale the main backgrounds by the factor corr discussed in the previous paragraph once it is hard to evaluate the impact of cuts on the other subleading backgrounds.
Beside the main backgrounds, the signals were rescaled by a conservative K-factor of 2 [58], irrespective the heavy scalar mass.

Discussion of the results
The missing transverse energy is the most effective variable to eliminate these backgrounds.
We scanned the E T cut from 200 up to 1500 GeV in order to investigate the best cut. The E T distribution of some representative signals and the main backgrounds and their cut efficiencies are shown in Fig. (5). We observe that the peak of the signal distributions scale roughly with the half the mass of the S scalar but the backgrounds are concentrated around E T ∼ 200 GeV.
For each point of the parameters space that passed all the constraints, we calculate the signal significance, N σ , with a metric that interpolates the gaussian and poissonian regimes at the same time it allows for including systematic effects in the background normalization without overestimating the significance [59,60] fixing a 10% systematic uncertainty in the total background normalization which seems to be realistic in view of the estimates from Ref. [54].
The missing E T cut was chosen in order to maximize the number of points of the parameters space with N σ ≥ 2, that is it, the maximum number of points that can be probed at the LHC at In the right panel, we show the cut efficiency in terms of the E T threshold. 95% CL. We chose to impose on the missing transverse energy of the events to get a more reliable estimate of the remaining backgrounds in the tail of the E T distribution. With this cut, around 26% of the parameters space that fulfills the constraint requirements have N σ ≥ 2.
The heat maps shown next display the density of points in selected bidimensional slices of the parameters space with the following extra requirement: N σ ≥ 2, and eventually with N σ < 2 or N σ ≥ 5. From these plots we want to understand which regions of the parameters space are likely to be probed at the LHC. In colder regions, the chance that a point of the parameters space can be probed at LHC is higher than in hotter regions. Drawing contour lines where N σ = 2 would define the regions that could be probed by the LHC searches for monophoton signals. We found more convenient to show the density of points that can be probed though. It is important to notice that encountering a point outside these 2σ regions means that it is not possible to reach this level of statistical significance with 3000 fb −1 but if a given point is found inside a certain 2D-slice of the parameters space it does not guarantee a signal, it is necessary to check if it falls inside all the other probable 2D-slices and then compute its statistical significance in this search channel. Looking at just one particular bidimensional region might be misleading. Yet, once its impossible to display the whole parameters space, we hope that showing 2D-slices of the most important parameters is a good guide for model building and further searches using the EFT approach.
Once we set c ae = c aµ , 78 possible bidimensional slices of the parameters space can be displayed, in principle, but just a subset of those parameters are actually important to our purposes.
Concerning the parameters related to the S scalar, c BB and c W W are less important than m S and c GG which control the production cross section and one of the main decay modes of the scalar S, and f S that controls its branching ratio to ALP pairs. Therefore, we do not consider the c BB and c W W dimensions. In the case of the parameters related to the ALP, we consider all of them but k BB and k W W , which are combined in the new parameter which is effectively involved in BR(a → γγ) and that can probed directly in many experiments de- Out of the 15 bidimensional slices that can be formed with the six parameters chosen, we display those whose interplay we found more important to be understood. Before discussing these plots, however, we show in Fig. (6) the points distributed in the BR(a → γγ) × BR(a → NN ) versus Λ space. These are two key variables to understand which models can be tested in colliders. In the upper panels, we show the distribution of points with N σ < 2 at left, that is it, those points that cannot be probed at ∼ 95% CL at the LHC, and, at right, the points for which the LHC will present sensitivity with N σ ≥ 2. We can identify two clearly distinct regions. Points with N σ < 2 have typically BR(a → γγ) × BR(a → NN ) < 10 −3 , while those with N σ ≥ 2 have BR(a → γγ) × BR(a → NN ) > 10 −3 and are more concentrated in the Λ < 20 TeV region. In the lower panels we show the prospects for an N σ ≥ 5 observation, for 3000 fb −1 and N σ ≥ 2 observation after 300 fb −1 at the left and right panel, respectively. The densities are qualitatively the same of the N σ ≥ 2 case but, of course, the number of points with N σ ≥ 5 is smaller than N σ ≥ 2. Overall, 26% of the points have N σ ≥ 2 for 3000 fb −1 , while 14% have N σ ≥ 5 for the same luminosity. We display the reach of the LHC to probe the EFT model for various luminosities in  at right. The mass of the S scalar is concentrated in the ∼ 500-800 GeV range, the ALP mass above ∼ 100 MeV and the product of branching ratios not far from 1%. Note that above ∼ 400 MeV, the signal in monopohotons fades away. The preferred value of the product of the branching ratios of the ALP into photons and sterile neutrinos, by its turn, can be explained in view of the ALP lifetime constraint depicted in Fig. (4). The larger the ALP mass, the easier is avoiding the lifetime constraint but it cannot be smaller than two neutrino masses nor large enough to open decays into pions and muons in order that the monophoton channel becomes effective. Second, in the lower panels of Fig. (8), we show the correlation between the branching ratio of S → aa and the product BR(a → γγ) × BR(a → NN ), at left, and the parameters f S at right. This time, we learn that reaching the desired statistical significance is possible in regions where the S decays almost all the time into ALP pairs and, naturally, that is possible mainly in the region with large Finally, we discuss the effective coupling ALP-photon, k γγ , that can be probed by several experiments as discussed in the Introduction. The model considered here presents an extended spectrum compared to the axion and ALP models where these particles interact just with the photon or the SM gauge bosons, however, apart from collider searches, we do not expect that prospects from other types of experiments will get considerably impacted. We therefore assume that bounds and prospects from any experiment probing just the ALP-photon coupling applies to our case. For collider searches, on the other hand, we should expect that the future projections of reach and sensitivity get shrunk with the addition of new decay channels that compete with the photon channel.
In the upper panels of Fig. (10), the heat maps of k γγ versus c aN and m N are shown at the left and right panels, respectively. ALP-photons couplings from 10 −6 to 10 −3 GeV −1 will be probed TeV and, again, k γγ in the interval 10 −6 -10 −3 GeV −1 . Yet, we can observe points with very large Λ when k γγ ∼ 10 −5 -10 −4 GeV −1 and even smaller. When Λ gets large, BR(S → aa) saturates to 100% but it does not affect the branching ratio of the ALPs, then this tail towards large Λ is the branching ratio of S → aa compensating for the drop in BR(a → γγ).
The lower right panel Fig. (10) is of particular importance. These are two parameters that use to be directly constrained in many experiments hunting for axions and ALPs. In this particular EFT model with an sterile neutrino, the bulk of points that can be probed at the LHC at 95% CL In order to compare the prospects of the HL-LHC to probe the ALP-sterile neutrino interaction with other experiments devised to investigate the ALP-photon coupling, we overlay our results, the red shaded areas, from the lower right plot of Fig. (10) into the m a × k γγ plane displaying several other prospects from other experiments in Figs. (11), (12). In these plots, we adapted the definition of our k γγ to match those of the works where we took the plots from when necessary.
Reminding that this region corresponds to a confidence level of 95% approximately.
We expect that in the presence of a new decay channel, like into sterile neutrinos, the branching ratio into photons gets smaller. The prospects from other experiments shown in Figs. (11), (12) then are expected to shrink actually. Nevertheless, we think it is very interesting to compare the reach of the LHC with those experiments to demonstrate the usefulness of LHC searches for this model.
Interestingly, the HL-LHC sensitivity to monophoton events from S → aa → γ jet + E T is complementary to the other experiments, covering a sizeable region of the m a × k γγ space not reached yet as we see in Fig. (11). At the left panel of Fig. (11), taken from Ref. [22], prospects from beam dump and the PrimEx recast [26] experiments, and regions previously covered by the LEP search for photons pairs (in yellow) and supernova constraints (in green) are shown. In all these experiments, the rate at which ALPs decay to photons is an important parameter and suppressing it softens the current bounds and represent a limiting factor for the future projections.
In the case where new decay modes are open to the ALP, just like sterile neutrinos or anything else, we should expect that the reach of these experiments get diminished. Searching for alternative ways to observe ALPs thus becomes important.
In the right panel of Fig. (11), taken from Ref. [15], we see that S → aa → γ jet + E T is complementary to the other search channels at colliders. Contrary to channels that rely on large decay rates into photons pairs, associate decay channels might probe regions where ALP-photon coupling is smaller in comparison. More importantly, we might be missing ALP signals by looking for signals that require large ALP couplings to photons or gluons.
In Fig. (12), modes and large ALP-photon couplings are expected to produce no signals. In particular, ALPphoton couplings down to ∼ 10 −6 GeV −1 can be probed if the ALP decays to sterile neutrinos.
At the right panel, we show the reach of the γ jet + E T from the LHC in comparison to the Belle and BaBar reaches in dedicated channels, but now for an ALP with k BB = 0. Again, the LHC can probe smaller couplings to photons if the ALP decays to an sterile neutrino. As a final remark, notice that supernova constraints are at the verge of excluding ALP-photon couplings from this model.

PARTICLES AND STERILE NEUTRINOS
In this section, we present an UV complete model for ALP couplings to SM particles and sterile neutrinos that is effectively described by the EFT Lagrangian of Eq. (2) in the low energy limitwhen the heavy new fields can be integrated out. These limits are taken from Ref. [30] and updated with the SN 1987 limit from Ref. [22]. The prospects of the experiments will shrink with the addition of the sterile neutrino decay channel but the LHC one obtained in this work. At the right panel, g aW is the ALP-coupling assuming k BB = 0 with the limits taken from Ref. [27]. The points represent benchmarks of the model presented in Section 5.

ALP couplings to the SM gauge bosons and the heavy scalar
Let us first consider a generic renormalizable model that extends the SM by the addition of a complex scalar singlet field Φ ∼ (1, 1, 0) and a pair of left and right chiral fermionic fields where the numbers in parentheses are the usual field transformation properties under the SM gauge groups SU(3) C , SU(2) L , and U(1) Y , respectively. We assume an approximate global chiral U(1) symmetry with the remaining fields transforming trivially. Such a symmetry is taken to be broken explicitly in the scalar potential, V (Φ), and also spontaneously by the vacuum expectation value of the complex scalar singlet, Φ = 0. The Lagrangian is in which L SM represents the SM Lagrangian and the Yukawa coupling constant y is taken to be real. The covariant derivative acting on Ψ in Eq. (23) is given by where the generators of SU(3) C and SU(2) L in the Ψ representation are normalized according to depending on the representation dimension of the field. For the fermionic field in the fundamental representation of SU(3) C and SU(2) L we have Ψ ∼ (3, 2, Y Ψ /2), and it is defined that k Ψ C (3) = k Ψ L (2) = 1 2 . The scalar potential for the singlet field is where for simplicity the parameters are taken as all real, with m, µ 2 Φ > 0 and λ 1Φ , λ 2Φ > 0. The first and the last terms in Eq. (25) are the only ones we consider to break explicitly the U(1) symmetry in Eq. (22). Note that the Lagrangian in Eq. The scalar singlet field, with its vacuum expectation value Φ = v Φ / √ 2, is decomposed as and we identify the heavy scalar and the ALP fields with S(x) and a(x), respectively. Thus, minimization of the potential leads to the following quadratic masses for these fields We adopt the view that it is technically natural to assume S , by the reason that the limit m 2 → 0, λ 2Φ → 0 leads to an increasing of the symmetries of the model [5]. In the present case the symmetries are augmented by the U(1) chiral symmetry in Eq. (22), which turns out to be exact at the classical level.
Within the consideration that the symmetry in Eq. (22) is exact and if Ψ L,R transforms nontrivialy under SU(3) C , i. e., d Ψ C = 1, then we would have in Eq. (23) a sort of KSVZ axion model [61,62], which presents a solution to the strong CP problem through the Peccei-Quinn mechanism [1,2]. More precisely, the original KSVZ model is defined with the extra fermions as Ψ L,R ∼ (3, 1, Y Ψ /2). Due to the anomalous feature of the U(1) symmetry a(x) turns out to be the axion field, getting a mass m a = 5.70(7) meV 10 9 GeV fa [63,64]. For the specific KSVZ type model in Eq. (23) we have f a = v Φ . We see that our previous phenomenological analysis cannot contemplate an axion model of this type because it is not possible to have the axion mass in the interval 1 MeV ≤ m a ≤ 0.26 GeV for v Φ ≥ 10 3 GeV, which is the typical mass scale of S according to Eq. (27). Thus, we use the model in Eq. (23) as prototype for the ALP interactions only, taking m a and the scale v Φ as independent free parameters in Eq. (27).
With the vacuum expectation value Φ , the relevant terms from Eq. (23) giving rise to the effective couplings of the ALP are in the Lagrangian where Ψ = Ψ L + Ψ R and m Ψ = y √ 2 v Φ . The couplings of a(x) to the gauge fields are obtained through one loop triangle diagrams shown in Figure 13. The coupling of a(x) to the abelian U(1) Y gauge field B µ (x) can be straightforward computed compound an interaction vertex of −i m Ψ v Φ aΨγ 5 Ψ and two interaction vertices of g Y Ψ 2 Ψ / BΨ to form the one-loop triangle diagram. The result is in which α Y ≡ g 2 /4π, and the product d Ψ C d Ψ L of the Ψ L,R representation dimensions of SU(3) C and SU(2) L counts the number of degrees of freedom with hypercharge Y Ψ running in the loop.
For the effective coupling of a(x) to the nonabelian SU(2) L gauge fields W i µ (x) we have in which α 2 ≡ g 2 /4π and d Ψ C k Ψ L is the product of the SU(3) C degrees of freedom -the dimension d Ψ C -in the SU(2) L representation of Ψ L,R by the normalization coefficient defined with the trace Tr[τ Ψ i τ Ψ j ] = k Ψ L δ ij in the computation of the fermion loop in the diagram. Similarly, the effective coupling of a(x) to the nonabelian SU(3) C gauge fields G a µ (x) is in which α s ≡ g 2 s /4π. Observe that a(x) does not couple to gluons if Ψ L,R is a color singlet once in that case k Ψ C = 0. For the scalar field S(x) its effective couplings to the gauge bosons at the lowest order are also obtained through computation of one loop diagrams similar to the ones shown in Figure 13. The effective interaction Lagrangians of S(x) with B µ (x), W i µ (x) and G a µ (x) are given by The Eqs. Details of the calculations of loop diagrams using both the linear and the polar representation of Φ for similar models can be found, for example, in Ref. [65].
We can now identify the coefficients of S(x) and a(x) couplings to the gauge bosons in the effective Lagrangian of Eq. (2) according to (with α Y ≈ 0.009, α 2 ≈ 0.03 and α s ≈ 0.1): The trilinear interaction Saa, allowing the decay S → a a, arises directly from the scalar potential. Taking into account Eq. (25) we have that the coupling is proportional to the vacuum expectation value v Φ , i. e., From this we see that in Eq. (2) the trilinear coupling is f S = (λ 1Φ − 3λ 2Φ )v Φ .

ALP couplings to the sterile neutrinos and the seesaw mechanism
We now discuss the light sterile neutrinos in which the ALP could decay into. Although this is a fine tuning, we note that it is not much worst of what occurs for the electron and the light quarks in the context of the SM.
It is possible to have sterile neutrinos with masses at the 1-10 2 MeV scale resulting from a mass generation mechanism leading to sub-eV masses for the active neutrinos. In order to see this we add, besides the fermionic field Ψ, two sets of sterile neutrinos fields: N kL ∼ (1, 1, 0) and S kL ∼ (1, 1, 0), with the index k = 4, 5, 6. Taking into account also the SM lepton multiplets, L l ∼ (1, 2, −1/2) and l R ∼ (1, 1, −1), with l = e, µ, τ , it is assumed the Z 4 symmetry presented above under which Along with Φ and Ψ, we also assume that the field N kL is charged under the approximate U(1) symmetry transforming as N kL → e −iα N kL . This symmetry is taken to be broken in the interaction terms involving the SM fields L l , H (the Higgs doublet) and the singlet N kL . Thus, the renormalizable Yukawa Lagrangian involving the new singlet fields and invariant under the Z 4 symmetry is where:H ≡ H * , with the antisymmetric matrix 12 = − 12 = 1; y kl and g k k are 3 × 3 complex matrices; and the 3×3 matrix M Sk k is taken diagonal with entries larger than v Φ . We have defined the right-handed field through the conjugation operation such that S c kR ≡ (S kL ) c . Under the assumption that M Sk k , v Φ > v W ≈ 246 GeV it is reasonable to consider a seesaw mechanism with two stages. The first stage involves N kL and S kL , with the mass matrix for these sterile neutrinos having the typical texture of the canonical seesaw mechanism which is, in the basis (N kL S kL ), where M Dk k = g k k √ 2 v Φ . Block diagonalization leads to the approximate states where B = M † D [M * S ] −1 , with N kL and S kL having the following mass matrices Taking v Φ ≈ 10 3 GeV and M S k k = n k k V /2 we have For g k k , n k k of order one we could have sterile neutrinos with masses 1-100 MeV for 10 9 GeV ≥ V ≥ 10 7 GeV. In the second stage the Higgs doublet gets a vacuum expectation value H = [0 v W ] T . The fields S kL are much heavier than N kL and the active neutrinos so that they can be integrated out leading to the following mass effective Lagrangian in which the new 3 × 3 matrices are defined according to The mass matrix involving the active, ν lL , and the sterile neutrinos, N kL , in the basis (ν L N L ) is Again, performing an approximate block diagonalization of M νN the flavor states of the active neutrinos ν lL mix with the sterile neutrinos states N kL to form the block diagonal states in which b = m † D M * −1 N is the mixing matrix; with N kL and ν lL having the following approximate mass matrices A discussion about the diagonalization procedure for the neutrinos mass matrices, including higher order corrections, involved in several seesaw mechanisms can be found in [66].
On the other hand, by the fact that the term y kl N c kRH † L l in Eq. (38) breaks the U(1) symmetry it is technically natural to have y kl 1. The mass scale of m L is then well below to the one required to be the compatible with the neutrinos mass square differences inferred from the experiments.
Therefore, we disregard m L so that the active neutrinos states ν lL have the mass matrix Depending on the scale V it would be possible to obtain active neutrinos mass eigenstates with mass at the sub-eV scale if y kl ≤ 10 −7.5 , corresponding to V ≥ 10 7 GeV. As said before, such small values would reflect the fact that the U(1) symmetry is almost exact.
The mixing matrix of the active and sterile light neutrinos has entries such that For the considered range of values for V this implies we estimate the active sterile neutrinos mixing for the particular model considered here to be such that |b kl | 2 ≤ 10 −5 (V ≤ 10 9 GeV).
This estimation can be compared with the limits from laboratory searches on sterile neutrinos through their mixing with active neutrinos. Denoting by N k the sterile neutrinos mass eigenstates of Eq. (47), with the lightest one being N 4 with mass in the range 1 MeV-100 MeV, the limits on the mixing with electron neutrino has upper bounds which varies from |U e4 | 2 4 × 10 −4 to |U e4 | 2 10 −8 at 90% C.L., as can be seen in Refs. [53,67]. These upper bounds could be compared to |b † e4 | 2 from the above assuming this is the mixing of the electron neutrino with the lightest sterile neutrino N 4 . For the mixing of the sterile neutrino, in the same mass range, with the muon neutrino the actual upper bounds varies from |U µ4 | 2 4 × 10 −2 to |U µ4 | 2 few × 10 −5 at 90% C.L. [53,67]. Thus, the model example we present here, with its mechanism for generating mass at the eV and MeV scales, can be compatible with the direct search experiments on sterile neutrinos.
There are cosmological constraints over sterile neutrinos that depend on their masses, timelife and mixing angle with the active neutrinos. Within the standard cosmological scenario, these constraints can be stronger than those from direct laboratory searches for sterile neutrinos that decay during or after the Big Bang Nucleosynthesis [53,68,69]. If sterile neutrinos produced in the primordial Universe have a lifetime well below 1 s they decay before the Big Bang Nucleosynthesis.
In this circumstance, their decay products might have reach thermal equilibrium with the other particles not affecting the Big Bang Nucleosynthesis. It is possible for to have lifetimes well below 1 s, as we estimate below. As far as we know, the actual constraints from cosmology do not consider the ALPs concomitantly with sterile neutrinos. Thus, an additional study is required to precisely define the extension of the cosmological constraints to a scenario like the one we take into account here.
The second term in Eq. (38) leads to the interactions where N k are the mass eigenstates of M N k k in Eq. (47). If the ALP can just decay into the lightest sterile neutrino N 4 ≡ N then we have from Eq. (2) an interaction of the same form of the effective Lagragian in Eq. (2) for describing the ALP-sterile neutrino interaction. From Eq. (52) we identify the effective coupling in Eq. (2) according to Observe that typical values c aN = 50 and Λ = 100 TeV give in fact v Φ = 1 TeV. Thus, with the potential in Eq. (25) and the Lagrangians in Eqs. (28) and (38) the effective Lagrangian in Eq.
(2) can be generated considering simple models which could be part of a complete theory at the ultra-violet regime.
It is important to observe that the lightest sterile neutrino N produced in the ALP decay, a → N + N , studied here might have a timelife long enough to decay outside the LHC detectors.
The main decay channels available for the sterile neutrinos with mass below that of the pion are: N → e − + W * + → e − + e + + ν; N → ν + Z * → ν + e − + e + ; and N → ν + Z * → ν + ν + ν. These From Eq. (36) we have f S ≈ m 2 S /2v Φ = 1 TeV. This model corresponds to the effective theory in Eq. (2) with c ae = c aµ = 0. Also, from Eqs. (21) and (55) we have which depends only on v Φ and the representation of the fermionic field Ψ. Some benchmarks are shown in Table (2) for some choices of representation of Ψ. As examples we consider: a quark singlet with electric charge 2/3, Ψ L,R ∼ (d Ψ C , d Ψ L , Y Ψ /2) = (3, 1, 2/3); a doublet of quarks, Ψ L,R ∼ (3, 2, 1/6); and a triplet of quarks, Ψ L,R ∼ (3, 3, −1/3). These benchmarks are represented as points in Figs. (11) and (12) in the ALP-photon coupling versus ALP mass plane. They are right in the region where the LHC is able to probe the EFT parametrization of the model. The signal significance of these points are above 5σ for 300 fb −1 . This is a consequence of having suppressed couplings to leptons that increases the product BR(a → γγ) × BR(a → N N ). Not only this, but an 1.4 TeV S scalar has a favorable cut efficiency compared to the backgrounds which are very efficiently eliminated by an 1 TeV cut in the missing transverse momentum.
It is worthwhile to observe that once the parameter f S scales with m 2 S 2λ 1Φ v 2 Φ , the decay width of S → a a behaves as Γ(S → aa) m 3 S 32πv 2 Φ and the branching ratio BR(S → aa) might increase with m S favoring the resonant ALP production.

CONCLUSIONS
ALPs are predicted in many extensions of the SM and their discovery could bring light to many problems in physics, astrophysics and cosmology. Sterile neutrinos, by their turn, might well be the first particles beyond the SM paradigm to be discovered if the signals from neutrino baseline experiments are confirmed. On the hand, an extended scalar sector with more scalars would be plausible to give mass to an extended matter sector. These three new ingredients thus sum up to a good target to phenomenological studies and model building.
In this work, we have investigated an effective field theory where a heavy scalar and an sterile neutrino couple to an ALP. We provided prospects to probe regions of the EFT parameters space at the 13 TeV LHC after 3000 fb −1 respecting several experimental constraints. The LHC is able to constrain and discover such models in monophoton searches with large missing transverse energy for ALP masses from 10 to 400 MeV, sterile neutrino masses from 100 keV to 100 MeV and S scalars from 500 GeV to 2 TeV. Effective ALP-sterile neutrinos and ALP-photons couplings, c aN /Λ and k γγ , as small as 10 −5 GeV −1 and 10 −6 GeV −1 , respectively, can be probed. Comparisons with existing constraints and prospects from other experiments show that the LHC in the monophoton channel can probe regions of the parameters space not accessible by those experiments establishing its complementary role in the quest for these kinds of physics beyond the SM as shown in Figs. (11) and (12).
To give a concrete example in model building, we present an ultra-violet complete model leading to the effective low energy Lagrangian studied in Section 2, where the ALP couples to sterile neutrinos whose masses can be well within the 1-100 MeV scale. The model contains a seesaw mechanism, which involves the scale of spontaneous breakdown of the approximate U(1) symmetry associated to the ALP, for generating mass for both sterile and active neutrinos. We found that singlet, doublet and triplet matter representations can be discovered after 300 fb −1 by the proposed strategy and respecting all the constraints considered in the EFT approach if the ALP couplings to charged leptons are suppressed.
Further research along the lines we have develop here could aim, for example, the sterile neutrino as a dark matter candidate, or a scenario with a few more sterile neutrinos with non-negligible mixing with the SM neutrinos generating couplings of type gN ν e and that, just like the models of Ref. [38,39], could explain excesses in the LSND and MiniBooNE data. Concerning this last proposed scenario, an UV complete model addressing the neutrino data could be built with the model presented in this work as an starting point.