A 750 GeV Portal: LHC Phenomenology and Dark Matter Candidates

We study the effective field theory obtained by extending the Standard Model field content with two singlets: a 750 GeV (pseudo-)scalar and a stable fermion. Accounting for collider productions initiated by both gluon and photon fusion, we investigate where the theory is consistent with both the LHC diphoton excess and bounds from Run 1. We analyze dark matter phenomenology in such regions, including relic density constraints as well as collider, direct, and indirect bounds. Scalar portal dark matter models are very close to limits from direct detection and mono-jet searches if gluon fusion dominates, and not constrained at all otherwise. Pseudo-scalar models are challenged by photon line limits and mono-jet searches in most of the parameter space.


Introduction
The ATLAS [1] and CMS collaborations [2] recently pointed out an excess in diphoton events with a peak in the invariant mass distribution around m γγ 750 GeV. Upon interpreting the events as the production and two-body decay of a new 750 GeV particle, current data cannot discriminate between a narrow or broad (up to ∼ 45 GeV) resonance. Although the evidence is far from conclusive, if it is confirmed with more luminosity it would be a monumental discovery after decades of undisputed success for the Standard Model (SM). Furthermore, it is natural to believe that such a hypothetical particle is linked to a bigger framework addressing, for istance, the gauge hierarchy problem and would be the herald of additional discoveries.
A more robust and elderly motivation for physics beyond the SM is the evidence for dark matter (DM) [3]. Among several candidates [4], a weakly interacting massing particle (WIMP) produced through thermal freeze-out [5][6][7] is undeniably one of the most appealing. Thus it is tempting to investigate whether the potentially new 750 GeV degree of freedom could act as a portal field, allowing DM and the SM to communicate beyond gravitational interactions.
This work focuses on (pseudo-)scalar portals and fermion DM candidates, both SM singlets. New (peudo-)scalars are ubiquitous in well-motivated frameworks for physics beyond the SM. At the same time, fermion singlets are DM candidates begging for new weak scale degrees of freedom, as gauge invariance forbids renormalizable interactions with SM particles [8]. We work within an Effective Field Theory (EFT) framework and write down the minimal theory for the LHC diphoton excess with a DM candidate. We define the theory at a cutoff scale Λ interpreted as the scale where heavy degrees of freedom are integrated out and we apply EFT methods to connect the interactions at different scales. While we present our analysis and results for the case of a Dirac fermion DM, it is straightforward to generalize to the Majorana case.
We start our phenomenological study with a comprehensive analysis of LHC results. Two different mechanisms, gluon and photon fusion, can be responsible for the (pseudo-)scalar production at colliders. In spite of being mediated by strong interactions, gluon fusion does not necessarily have to be the dominant production mechanism at the LHC since we have no actual evidence that the new particle couples to gluons at all. From the diphoton excess, we do know that the resonance must couple to photons. This implies that there exists an irreducible photon-fusion contribution to the resonance production, which can be dominant one or not depending on the relative sizes of the couplings to photons and gluons. We therefore include both production mechanisms in our study, and we identify where the EFT is capable of accounting for the diphoton events while at the same time being consistent with √ s = 8 TeV data. The presence of a DM candidate in the EFT impacts our analysis even before discussing any DM phenomenology. Once produced at the LHC, the (pseudo-)scalar could be allowed to decay to invisible final states, altering the width and diphoton rate. For this reason, we find it convenient to divide our LHC study into two scenarios: SM Dominated Resonance: The DM mass is above the critical value 375 GeV. The resonance only decays to SM final states and it is typically narrow. The ATLAS preferred value of Γ S 45 GeV can be obtained only for large couplings to SM fields which are inconsistent with searches in other decay channels such as Zγ.
DM Dominated Resonance: The DM mass is below 375 GeV such that decays to DM pairs are kinematically allowed. This invisible channel is very likely to dominate the total width and the resonance is now quite broad.
In each of the above we perform a thorough exploration of the parameter space. The presence of a sizeable coupling to gluons utterly drives LHC phenomenology, as gluon fusion is clearly the leading candidate for the (pseudo-)scalar production. However, photon fusion can still dominate the production if the coupling to gluons is small enough. For example, this would be the case for UV-complete theories where any heavy particles that are integrated out at the cutoff scale Λ do not carry color charge. Despite the apparently large parameter space, we identify two main EFT regimes where the production is dominated by a single partonic process and where the couplings of the new particle to SM gauge bosons are quite constrained. We emphasize that every specific UV completion with no additional degrees of freedom below the cutoff Λ must satisfy the constraints of our EFT analysis.
In the second part of our study we incorporate the DM phenomenology. For parameters favored by LHC data, we further impose constraints from DM searches and also identify regions where the DM has a correct thermal relic density. Collider searches for mono-jet events turn out to be relevant only in the DM dominated scenario, as the associated cross section falls rapidly as the DM mass increases above the resonant value 375 GeV. Direct Detection (DD) experiments constrain only the scalar portal case, and the coupling to gluons is again crucial. If such a coupling is present, DD rates are dominated by the gluon content of the nucleons. If not, both the coupling to photons and the loop-induced couplings to gluons and light quarks contribute to the signal. In each case we evaluate DD rates through a rigorous Renormalization Group (RG) procedure, which is mandatory as the scale separation between DD and LHC experiments is large. On the contrary, indirect detection (ID) experiments put bounds only on pseudo-scalar mediators because annihilations mediated by a scalar portal are p-wave suppressed. This difference also explains why larger couplings to the scalar are necessary to reproduce the observed DM abundance through thermal freeze-out.
The paper is structured as follows. In Section 2 we introduce the EFT that will be the basis of our study. Section 3 deals with the connection between energy scales. In Section 4 we introduce the two different LHC scenarios and identify the parameter space region allowed by collider searches in both cases. Finally, we present our DM analysis in Section 5 and summarize our main findings in Section 6. We provide appendices with explicit expressions for decay rates and cross sections, details of the RG procedure, and methods for the relic density calculation.

The EFT for Dark Matter and the Diphoton Excess
We introduce the minimal EFT necessary to describe the diphoton excess at the LHC, while simultaneously providing a stable DM candidate. We augment the SM by two singlet fields: a real scalar S with mass m S = 750 GeV and a fermion χ. The formalism developed in this Section is valid for both Dirac and Majorana χ. Although we give the details of the EFT for the case of a scalar S, the generalization to the case of a pseudo-scalar P is straightforward as shown at the end of this Section.
Within our framework, the LHC excess is accounted for by the production of S and its subsequent decay to photons. At the same time, S also acts as a portal to the DM particle χ assumed to be a stable field as a consequence of a Z 2 symmetry. The EFT Lagrangian reads with L SM the SM Lagrangian. We organize the interactions in L int by distinguishing between renormalizable and non-renormalizable operators, and we further classify the latter according to their mass dimension The sum on the right-hand side of the above equation runs over all SM gauge invariant operators for each mass dimension d. The higher-dimensional operators O (d) α d are suppressed by powers of the EFT cutoff Λ, understood as the mass scale where heavy degrees of freedom generating the interactions are integrated out. The dimensionless and renormalization-scale, µ, dependent Wilson coefficients, c α d , encode unresolved dynamics above the EFT cutoff.
The renormalizable piece contains the portal interaction between the two singlets Additional renormalizable interactions in the scalar potential are not forbidden by any symmetry. In particular, operators involving both S and the SM Higgs doublet H would induce a mixing between S and the SM Higgs boson h, affecting production and decays of both. This scenario is quite constrained by Higgs coupling measurements and it has been recently studied in Refs. [9,10]. In this work, we assume these scalar potential interactions to be absent, as realized in several UV completions (see e.g. Ref. [11]). Moving on to higher-dimensional operators, we consider the d = 5 contact interactions Here, we assume the EFT cutoff to be much higher than the weak scale, Λ m Z , and therefore we couple our degrees of freedom in a SU (3) c × SU (2) L × U (1) Y gauge invariant way. Also in this case, other operators in Eq. (4) are in principle allowed by symmetry considerations: the Higgs portal operator χχH † H, the coupling to SM fermions SHf L f R , and additional d = 5 scalar potential interactions. We assume again that these couplings are not present at the EFT cutoff, as it would be the case in several UV completions (also discussed in Ref. [11]). However, assuming that they vanish at the cutoff does not save us from having them at other scales: the absence of a symmetry protection allows the RG evolution to switch them on through radiative corrections, and it the next Section we quantify how this happens. These radiative contributions play no role for LHC physics, and we can safely use the Lagrangian in Eq. (4) to study LHC phenomenology. The situation is rather different for DD, since we evolve the EFT all the way down to the scale of nuclear physics.
We stop at dimension 5 and we do not consider other operators in the sum in Eq.
(2), consistently with the fact that the EFT cutoff is taken to be much larger than the typical energy scale for LHC processes. To summarize, the EFT obtained by adding S and χ to the SM and with the interactions in Eq. (3) and Eq. (4) will be the basis for this work.
The EFT where the new bosonic degree of freedom is a pseudo-scalar P is very similar Unless we have a substantial amount of CP violation, we have either the scalar S or the pseudoscalar P . RG effects turn out to be negligible for the pseudo-scalar because DD constraints are very weak, so the results in Section 3 are relevant for the scalar only. Furthermore, LHC rates are identical for the two cases, so the analysis performed in Section 4 is valid for both. The DM phenomenology is drastically different between the two cases, since DM annihilations mediated by a (pseudo-)scalar are (s-)p-wave processes, and DD constraints are negligible for the pseudo-scalar. For this reason, we keep the DM discussion in Section 5 separated for the two cases.

RGE Scale Connection and Direct Detection Rates
The Wilson coefficients in Eq. (4) are generated at the cutoff Λ by integrating out heavy degrees of freedom, while LHC data bound their values at the typical collider scale. In order to perform a consistent EFT analysis, we would have to RG evolve the interactions down to LHC scales before putting limits. As we will show shortly, these corrections turn out to be inconsequential. Nevertheless, our EFT looks very different at energy scale of the order of ∼ 1 GeV, where nuclear matrix elements are evaluated to compute DD rates. This procedure can significantly affect DD rates, as pointed out for several cases in Refs. [12][13][14][15][16][17][18][19][20][21][22][23].
The only SM fields accessible at such a low-energy scale and relevant for DD observables are light quarks, gluons, and photons. We define a different EFT for DD in terms of these light degrees of freedom, and the relevant interactions for our study are the following with Wilson coefficients evaluated at the nuclear scale µ N 1 GeV. Our goal here is to connect the Wilson coefficients at the cutoff scale appearing in Eqs. (3) and (4) with the ones at the nuclear scale in Eq. (7). This is achieved by performing the RG evolution (RGE) obtained via the following steps • perform the RGE from µ = Λ down to the scalar mass µ = m S ; • integrate out the scalar field S at the scale µ = m S ; • perform the RGE from µ = m S down to the weak scale µ m Z ; • integrate out the heavy SM degrees of freedom (top, W, Z, Higgs); • perform the RGE from µ = m Z down to the nuclear scale µ µ N , and in the process integrate out the intermediate heavy quark (bottom and charm) at their mass threshold.
Here, we present the main RGE results with details of calculations deferred to Appendix B.

Running from Λ to m S ?
The RG evolution to lower scales has two main effects: multiplying couplings by overall constants (self-renormalization) and inducing new interactions (operator mixing). Here, we inspect if the latter is ever relevant for LHC physics, as it is the only process which could induce a different phenomenology. Self-renormalization can always be taken care of by considering the Wilson coefficients at µ = m S .
If the scalar couples to gluons, QCD running induces couplings to quarks at the scale m S . This can be phenomenologically relevant only for the top quark, since the effect is proportional to the Yukawa coupling. If sizeable, this coupling can contribute to the total width of the scalar and open thett production channel at the LHC. This effect is quantified by the induced partial width tott in units of the one to gluons where the Wilson coefficient ct t (m S ) can be obtained from the results in Appendix B This effect is too small to play any role at the LHC. The operator mixing and the radiatively induced interactions for the case of no coupling to gluons are even more suppressed as a consequence of the weak fine structure constant and smaller anomalous dimensions. Other potentially relevant RG effects arise from inducing operators involving both the new resonance and the the SM Higgs doublet. In our case we would induce the dimension 5 operator S(H † H) 2 , which has two main effects. First, it opens up the new S decay channel into two or more Higgs bosons. Second, it induces a mixing between S and h with consequent change of production and decay rates for both scalars. In particular, the couplings between h and other SM fields would be different with respect to their SM value. The Wilson coefficient for this dimension 5 operator as induced at the scale m S can be calculated using the analysis in Appendix B, and it results in Contributions from the gluonic coupling c GG only appear at the two-loop level and we neglect them. The decay width into Higgs bosons, relative to the one into photons, is then given by where c γγ (m S ) s 2 w c W W (Λ) + c 2 w c BB (Λ) (see Appendix A). This RG induced decay width is too small to play any role in our analysis. The induce mixing between S and h is quantified by the mixing angle with c H (m S ) still given in the expression in Eq. (11). As discussed extensively in Section 4, the coupling of the new resonance to electroweak gauge bosons can be at most c BB /Λ c W W /Λ = 0.3 TeV −1 in the photon fusion regime. Such couplings give rise to a very small mixing angle, α 10 −4 , well below any experimental constraints [9].

Connecting m S to µ N I: Scalar coupled to gluons
If c GG does not vanish at the cutoff, the couplings to electroweak gauge bosons provides a negligible contribution to DD rates, thus we ignore their effects and only consider QCD running. This has been extensively studied in the literature (see e.g. Refs. [17,22]). We repeat the leading order (LO) analysis for completeness in Appendix B, where we argue that next-to-LO corrections only modify the final result by a few percent. We summarize here the main results. As discussed in Section 3.1, we should only be concerned about the RG from m S to µ N . Thus we start our analysis at the scale m S , where we integrate out the scalar and write down the effective Lagrangian valid for energy scales above the top mass. The coupling is obtained via a tree-level matching The connection between the Wilson coefficient in Eq. (14) evaluated at the renormalization scale µ = m S and the ones of the effective Lagrangian for DD in Eq. (7) evaluated at the nuclear scale is achieved as follows 3.3 Connecting m S to µ N II: Scalar coupled EW gauge bosons On the other hand, if S does not couple to gluons at the scale Λ the running driven by electroweak gauge bosons turns out to be relevant for the rate calculation [14,20]. We start at the scale µ = m S , where we integrate out the scalar and end up with the effective Lagrangian The tree-level matching in this case is analogous The connection between the couplings in Eq. (18) and the ones for DD in Eq. (7) reads

RGE Analysis: Summary
If it useful to summarize the main results of this Section. The RGE from Λ to m S does not affect the LHC phenomenology, thus we use the EFT defined at the scale µ = m S to perform the LHC phenomenological analysis. Then we have to connect the couplings at µ = m S with DD rates. For coupling to gluons we use the results in Eqs. (16)- (17), whereas for interactions with electroweak gauge bosons we have the low-energy couplings given in Eqs.

Two Different Scenarios for LHC
The Lagrangians introduced in Section 2 contain seven free parameters: three mass scales (m S , m χ , Λ) and four dimensionless couplings (c GG , c W W , c BB , c χS ). We set m S = 750 GeV motivated by the diphoton excess. As justified in Section 3.1, we start with the EFT defined at µ = m S and present the results of our LHC analysis in terms of c XX /Λ, where X = {G, W, B}.
It is worth pointing out that we always have in mind values Λ few TeV to safely satisfy the EFT hypothesis. DD rates are computed through the RGE from m S down to the nuclear scale as discussed in Section 3, thus they also depend on the combination c XX /Λ only. This leaves us with the DM mass and four couplings. The DM mass value m S /2 375 GeV is quite special, as for masses smaller (greater) than this critical value the scalar S is (not) allowed to decay to DM pairs. In view of this, we divide our study into these two main cases. They correspond to quite different scenarios at the LHC, since the scalar resonance is typically narrow unless we open the decay to DM. The origin of this can be traced back to the fact that decays to DM are the only ones mediated by a renormalizable interaction. Before introducing and studying the two scenarios, we discuss the cross sections for gluon and photon fusion. Related LHC studies considering gluon fusion partonic processes have been recently performed in Refs. [9,10,[24][25][26][27][28][29][30][31][32][33][34][35][36][37][38], whereas photon fusion was considered in Refs. [38][39][40][41][42]. If the coupling c GG is non vanishing at the cutoff then gluon fusion certainly dominates partonic productions for every channel. In the absence of such a coupling at the cutoff scale, one may wonder about the main production mechanism. The two-loop induced coupling c GG at the scale m S turns out to induce a gluon fusion rate that it is subdominant compared to the vector boson fusion (VBF) contribution. Strictly speaking, all possible VBFs contribute to the cross section in our EFT, namely partonic processes with initial state ZZ, W W , W Z, Zγ, W γ and γγ. In particular, we cannot have only the γγ process since this would require to have the three effective vertices vanishing (SW W , SZZ, and SZγ) with the only freedom of tuning the two Wilson coefficients c BB and c W W , as pointed out in Ref. [38]. Photon fusion diagrams dominate at LO, and the next relevant contribution is the interference between photon and weak boson processes. We neglect this correction since it is approximately only a 10% modification of the total cross section [40]. However, these couplings to weak gauge bosons lead to proton-proton collisions with W W , ZZ, and Zγ final states that are bounded from LHC Run 1 searches.
We include both gluon and photon fusion in our analysis and identify for what parameters each one is dominant. As already stressed in Section 2, the LHC analysis performed here is valid for both scalar and pseudo-scalar mediators.

LHC production cross sections
Whether the resonance is broad or not, the question about the partonic production mechanism is still open. The general formalism for LHC cross sections can be found in Appendix A. Here, we specialize the general expression given in Eq. (87) to the two cases of our interest.
Gluon fusion is a natural candidate, as long as the coupling to gluons c GG is switched on. If this is the case, the production cross section results in The overall multiplicative factor K GG accounts for QCD higher-order contributions. This Kfactor correction does not depend on the CM energy of the proton-proton collision, as we are always interested in resonant processes, and in what follows we use the value K GG = 1.48 [26]. The quantity C √ s GG is defined in Eq. (88) as an integral over the gluon parton distribution function (pdf) in the proton. We adopt the gluon pdf from Ref. [43], and using their public code 1 we find the following numbers Thus cross sections at 13 TeV are rescaled from the ones at 8 TeV by the following factor Partonic productions through photon fusion have a cross section The factor r inel accounts for inelastic processes where the proton gets destroyed after radiating a photon. Unfortunately, its precise value suffers from theoretical uncertainties. The recent LO calculation with Madgraph [44] performed in Ref. [41] found that the elastic processes are only 4% of the total events, or equivalently r inel = 25. This is consistent with the discussion in Ref. [40], claiming the range r inel ∈ [15,25]. We normalize our √ s = 13 TeV cross section with the results of Ref. [41]. Upon setting r inel = 25, we find the contribution from elastic processes A key quantity in this regime is the rescaling between cross sections The output of the Madgraph calculation in Ref. [41] gives R γγ 2. This quantity, however, is not very well known for instance due to uncertainties regarding the inverse proton radius and the size of r inel . Following the discussions in Refs. [39][40][41], we take the range R γγ ∈ [2, 5], and present our results for the two representative values 2 The expressions in Eqs. (25) and (29) are the master equations for this Section. Combined with the decay widths listed in Appendix A, they allow us to compute the production cross section for any ij-pair. We use them to find the EFT parameters consistent with both the diphoton excess [1, 2] and the LHC Run 1 constraints listed in Tab. 1. For the signal we identify the 1σ and 2σ regions consistent with the cross section We do not consider limits fromtt searches [46,47] since the rate is loop-suppressed in our EFT.
In addition, we do not show in our plots γγ limits at √ s = 8 TeV. They are definitely consistent with the diphoton signal at √ s = 13 TeV for the case of gluon fusion, given the rescaling factor in Eq. (28). The photon fusion presents tension for R γγ = 2, the lowest rescaling factor we consider, but is consistent for the larger value R γγ = 5. After this exploration of the parameter space, we will consider constraints from relic density and DM searches.

A SM Dominated Resonance
Our first scenario features DM masses above the kinematical threshold for decays m χ m S /2. The total width of the scalar S results in obtained by using the decay widths given in Appendix A in the m W,Z m S limit. In this scenario, the total width Γ S is quite narrow for a large region of the parameter space because of the (m S /Λ) 2 suppression, consequence of the non-renormalizable interactions in Eq. (4).
Gluon and photon fusions are both a potential source for production and which one dominates depends on the relative sizes of the couplings. As discussed in Section 4.1, the photon fusion cross section suffers from theoretical uncertainties in the inverse proton radius. We summarize our results in Fig. 1 and Fig. 2, where we choose R γγ = 2 and R γγ = 5, respectively. These plots show the same quantities with the only difference being the choice for R γγ , so we discuss each panel only once and we emphasize whenever the choice for R γγ makes any difference. Figure 1: Parameters for diphoton excess (green regions) and excluded by LHC Run 1 searches (red and blue regions). The Γ S ≥ m S region is completely shaded away, whereas the one with Γ S /m S ≥ 6% is shaded with light gray. We set the rescaling factor defined in Eq. (31) to R γγ = 2. In the upper panels we switch on the coupling to gluons and consider c W W = 0 (left) and c BB = 0 (right). In the lower panels we assume the production dominated by photon (left) or gluon (right) fusion and visualize the parameter space in the (c W W , c BB ) plane.
In both figures, the green shaded areas corresponds to the 1σ and 2σ regions for the diphoton excess. We also shade the areas excluded by LHC Run 1 searches with red (W W , ZZ, Zγ) and blue (di-jet) . Finally, we ignore parameters giving Γ S ≥ m S and shade with light gray the regions where Γ S /m S is above the 6% value favored by ATLAS. We start our discussion from the two upper panels, where in the left (right) we show the parameter space for only c GG and c BB (c W W ) switched on. At very small values of c GG , located on the left of the plots, the production is dominated by photon fusion where s w (c w ) is the sine (cosine) of the weak mixing angle. Not surprisingly, both the region accounting for the diphoton excess and the exclusion limits show up at constant values of c BB (left) or c W W (right). It is always possible to have a consistent explanation of the diphoton excess through photon fusion with only c BB , although the R γγ = 2 case features some tension with results from Zγ searches. However, the case with only c W W is excluded by Run 1 results. As we move toward the right of the two upper panels, eventually we increase c GG enough such that gluon fusion is the dominant production process while still being consistent with dijet searches. In this opposite limit the cross section is approximately The diphoton excess is again accounted for by a constant value of c BB (left) or c W W (right). Again, having c W W only is excluded by Zγ limits, whereas the case with only c BB is allowed regardless of the specific value of R γγ which plays no role for gluon-fusion dominated processes. A thorough exploration of these two opposite limits is provided in the lower two panels of Fig. 1 and Fig. 2. More specifically, we consider on the left (right) values of c GG small (big) enough such that the production is dominated by photon (gluon) fusion, and we visualize the allowed regions in the (c W W , c BB ) plane. First, we notice consistency with our previous findings. Regions with c BB c W W can account for the diphoton events, with again some tension for the case with a rescaling factor R γγ = 2. Conversely, the case with mostly c W W is badly excluded by Run 1 searches. An interesting intermediate case, allowed for both gluon and photon fusion, is for couplings to electroweak gauge bosons roughly of the same size c BB c W W . In particular, in the photon fusion case R γγ = 5, right at the edge of the Zγ limit the couplings c BB and c W W can be large enough to give a relatively broad resonance of Γ S /m S (2 − 3)%. However, it is not possible to reproduce the ATLAS preferred value Γ S /m S 6% as the Zγ limits are too stringent. Nevertheless, this seems to be the only point in parameter space where a sizable width is still possible without having decays into invisible states (as discussed below).

A DM Dominated Resonance
For DM mass values below m S /2 the LHC phenomenology is drastically different. The decay channel to DM is now open and completely dominates the total width The above equation is obtained by using the results of Appendix A and ignoring the phase space suppression, only relevant for DM masses very close to m S /2 350 GeV. The ATLAS best fit value for the witdh Γ S 45 GeV is easily obtained for couplings to DM of order one. The analysis proceeds similarly to the SM dominated resonance scenario in Section 4.2, with the important exception that we have also the coupling c χS in the game. As a consequence, we can always fix the total width Γ S to any value. Our figures in this Section follow the same conventions as the ones adopted in Section 4.2, with two important differences. First, we present our results in this DM dominated resonance scenario for the ATLAS best fit value Γ S /m S 6%, since we have the freedom to independently choose Γ S . We use arrows in our plots to show how our results change if one choses a different value (note that the limits from Zγ, ZZ, and W W searches also follow the arrows). Second, we shade with light gray in each plot the region below the center of green bands in Figs. 1 and 2. At the boundary of this "SM dominated" portion of the parameter space, the SM contribution alone accounts for the signal, therefore we cannot go below it. Figure 3: Parameter space for the DM dominated scenario for the rescaling factor defined in Eq. (31) equal to R γγ = 2. We identify the regions preferred by the diphoton excess (green) and excluded by LHC Run 1 (red and blue). We always set Γ/m S 6%, and the little arrows show how our bands moves as we change this value. We shade with dark gray the Γ S ≥ m S region, and with light gray the region below the boundary where the decay width to SM states alone accounts for the signal. In the upper panels we consider couplings to gluons and c W W = 0 (left) or c BB = 0 (right). In the lower panels we assume the production dominated by photon (left) or gluon (right) fusion and visualize the parameter space in the (c W W , c BB ) plane.
The results are shows in Figs. 3 and 4, where the only difference between the two figures is still the choice of R γγ . As usual, we start from the class of models where we only have the scalar The case of only couplings to c BB works in this regime, with again some tension with Run 1 bounds for R γγ = 2. The case with c W W only, other than being excluded (in agreement with Ref. [40]), also falls well inside the SM dominated region and therefore we neglect it. On the opposite end of the plots, gluon fusion dominates all productions with cross section Unlike the previous scenario, this gluon fusion regime does not pinpoint a specific value of c BB or c W W , but the green bands roughly corresponds to c BB c GG const, with this behavior persisting for values of Γ S /m S not too close to the SM dominated gray region. This is of course due to the fact that the total width is dominated by invisible decays, and not by decays to gluons as in the SM dominated case. The gluon fusion regime is again consistent with data if we only have the coupling c BB , whereas the case with only c W W is excluded by LHC bounds.
As done before, in the bottom panels of Figs. 3 and 4 we further explore the (c W W , c BB ) plane for the opposite photon (left) and gluon (right) fusion regimes. Not surprisingly, the allowed values are c W W c BB up to c W W c BB .

Dark Matter with a (Pseudo-)Scalar Portal
With the LHC analysis performed in Section 4, we are ready to include the DM in our study.
Recent and related DM works on the possibility of a 750 GeV (pseudo-)scalar mediator can be found in Refs. [11,[56][57][58][59][60][61][62][63][64]. 3 We extend their DM analysis by using the output of our LHC study, where we have considered the full parameter space with both gluon and photon fusion active. We present our results for both cases of scalar and pseudo-scalar mediator. It is useful again to divide our discussion between the two scenarios of SM and DM dominated resonance, presented in Sections 4.2 and 4.3, respectively. For each case we compute the relic density of the DM as a function of its mass by following the procedure outlined in Appendix C and we demand that it makes all of the measured DM in the Universe (Ω χ h 2 = 0.1188 ± 0.0010 as inferred by the latest results of Planck [71]). Furthermore we impose constraints from the following DM searches: • Collider searches for events with a singlet jet and missing energy (j+MET) [72,73], which are suitable only in the gluon fusion regime. We implement our EFT in FeynRules [74] and we generate the associated UFO model file [75]. The signal in then obtained by using MadGraph [76]. We impose the bound in the MET > 500 GeV bin, where the signal must satisfy the bound σ(j + MET) 6 fb. The DM analysis is performed by using the full results of our simulations. Here, we give the j+MET production cross section for two opposite limits in order to understand the qualitative behavior of this constraint. At DM masses well below the resonant value m S,P /2, the signal cross section depends only on the coupling to gluons. For CM energy √ s = 8 TeV it scales as 4 The mediator is produced on-shell in this regime and then it decays to DM pairs with 100% branching ratio. This explain the absence of c χS in Eq. (40), which holds as long as the DM coupling c χS is such that the scalar decay width is not too large (see Eq. (37)). We checked that Eq. (40) correctly describes the parameter space region we are interested in. Conversely, the mediator is way off-shell for DM masses above the resonant value m S,P /2. The process in this case can be approximately described by a contact interaction between gluons and DM particles [77]. For such a heavy DM particle we have As a consequence, collider limits do not play any role for the SM dominated scenario. The reader interested in further details can find a specific mono-jet analysis for the 750 GeV portal in Ref. [78].
• Direct searches, where we impose the most recent LUX bounds [79] and show LZ projections as extracted from Ref. [80]. These limits are only relevant for scalar mediators. Here, we present the spin-independent cross section for a DM Dirac fermion scattering elastically off a nucleus derived from the interactions in Eq. (7). These low-energy couplings are connected to those at the LHC scale as explained in Section 3. For the DM-quark and DM-gluon operators we follow the steps in Ref. [81], while for the DM-photon interactions we follow Ref. [82], and write where m A is the target nuclear mass. In the rest of the text we denote the scattering cross section of a single nucleon by σ SI dividing Eq. (42) for (A 2 F 2 H ) and replacing the DM-nucleus reduced mass with the DM-nucleon one. Here we have defined denoting, respectively, the contributions from C q , C G , and C F to the scattering amplitude. For the up and down scalar couplings we use the recent determination in Refs. [83,84] based on chiral perturbation theory and a Roy-Steiner analysis in good agreement with Ref. [85]. For the strange scalar coupling we use the lattice QCD calculation f p s = f n s = 0.043 ± 0.011 [86]. In the analysis below we use the central values for these matrix elements. These values then give the gluon coupling f N G = 0.894. b(A) is the harmonic oscillator parameter defined in Ref. [87] b(A) = 41.467 Finally, with F i F j we denote the nuclear form factor averaged over the velocity integral and the detector efficiency [82]. We follow the analysis of the LUX experiment and use a standard Maxwellian velocity distribution with v 0 = 220 km/s. The Helm (F H ) and Raleigh (F R ) form factors we take, respectively, from Refs. [87] and [82] where, for the latter, we set the two-body parameter c 2 = 0.
• Indirect searches, which on the contrary are only relevant for pseudo-scalar mediators since DM annihilations mediated by a scalar field are p-wave suppressed. We impose limits on the annihilation cross section from γ-ray line searches from both Fermi [88] and H.E.S.S. [89] considering peaked and cored DM density profiles, as well as limits on dwarf spheroidal galaxies (dSphs) observations by Fermi [90]. Notice that the H.E.S.S. collaboration imposes limits in [89] only for an Einasto profile. Since the bound from H.E.S.S. are very sensitive to the choice of the profile (the region of interest is a small circle of 1 • centered in the galactic center), we also consider a cored profile (Burkert) in our analysis by rescaling their limits with the J-factor given in Ref. [92]. We compare the experimental bounds on the annihilation in lines (Fig. 8 of [88] and Fig. 4 of [89]) with the predicted annihilation cross section into γγ + γZ/2. Furthermore, when imposing continuum limits from dSphs we take advantage of the two following facts [92]: photon spectra from electroweak gauge boson radiation is almost universal (in this case we compare the experimental bounds on W + W − (bottom right-panel of Fig. 8 of [90]) with the predicted annihilation cross section into W + W − + ZZ + γZ/2), as it is the case for the ones initiated by gluons and light quarks (in this case we compare the limits onūu (top right-panel of Fig. 8 of [90]) with the predicted annihilation cross section into gluons). We rescale all indirect limits by a factor of 2 to account for our choice of Dirac DM.

SM Dominated Resonance
We start by considering DM masses above the critical value m S /2, therefore the (pseudo-)scalar mediator can only decay to SM final states. A thorough exploration of the parameter space in this scenario was performed in Section 4.2, with regions consistent with LHC results shown in Figs. 1 and 2. We study DM phenomenology for three representative classes of models. We start by examining UV completions yielding only the coupling c BB to the hypercharge gauge boson. Our results are shown in Fig. 5. In the scalar mediator case (left), current and projected direct searches are not capable of probing the thermal relic region. In fact, they cannot probe any point of the region where the coupling c χS is perturbative, as the radiatively induced couplings to quarks and gluons given in Eqs. (21)-(24) are too small. DM production through thermal freeze-out can be analyzed in two different regimes. DM particles with mass in the range m S /2 < m χ < m S can only annihilate into SM fields. For DM masses away from the resonance this requires rather large couplings to the scalar portal, almost up to c χS = 4π for m χ m S , as a consequence of the p-wave suppression. Annihilations to mediators through the process χχ → SS open up for DM mass values above m S . The required value for c χS suddenly drops above this threshold, and it increases again for larger DM masses. We mention the tantalizing possibility of probing the scalar portal in the photon fusion regime through ID via the process χχ → SS → 4γ. In this case, the photons are distributed in a box centered around m S , and for DM masses not too much larger than m S they exhibit spectral features similar to the case of a line (see e.g. [91] for a dedicated study of γ-ray boxes with the forthcoming Cherenkov Telescope Array). We leave this direction for future work.
Results for the the pseudo-scalar case are shown in the right panel of Fig. 5. We observe a similar feature in the relic density line, although far less pronunced. The drop in c χP is much smaller because annihilation into SM fields is an s-wave process. This also implies that ID limits are very stringent. In the lower DM region (m χ 500 GeV) Fermi rules out thermal relics, whereas the H.E.S.S. line limits are excluding the thermal region only for the choice of an Einasto density profile.
A potentially interesting intermediate case is photon fusion at the LHC but with both couplings c BB and c W W present. As shown in Figs. 1 and 2, the most we can push is for c W W c BB with a small region where c W W can be a few times larger than c BB such that the scalar width is relatively broad. The conclusions for dark matter phenomenology are pretty similar to the case in Fig. 5. The only differences are that the relic line will move towards lower values of c χS and c χP . For the pseudo-scalar case, limits from line can be softened as a consequence of the continuum γ-ray contamination from c W W . These cases only form a small part of the allowed parameter space (see Figs. 1 and 2), and since we cannot push c W W much above c BB without running into conflicts with Zγ limits, we do not further discuss this case.
We now consider models where the (pseudo-)scalar couples to gluons and LHC productions are dominated by gluon fusion. Unlike the previous case, the couplings are not univocally determined, as gluon fusion dominates and di-jet constraints are satisfied in the whole range 0.03 TeV −1 c GG /Λ 0.07 TeV −1 . We choose the value at the center of this allowed range and we show results for this scenario in Fig. 6. The couplings to gluons for a scalar mediator is responsible for quite large direct detection rates. As an example, the RG analysis in Section (3.2) yields a direct detection cross section σ SI c 2 χS × 2.2 · 10 −46 cm 2 for m χ = 1 TeV and c GG /Λ 0.03 TeV −1 . Limits from mono-jet events are not relevant in this DM mass range (see Eq. (40)). The thermal relic line for m χ < m S is almost completely excluded by LUX, except for DM masses extremely close to m S /2. Similar to the photon fusion case, for m χ > m S the required value of c χS suddenly drops and a thermal relic is consistent with LUX bounds. However, the entire parameter space will be deeply probed by LZ. Although the results in Fig. 6 are presented for a fixed value of c GG /Λ, it is straightforward to rescale the results. DD bounds scale linearly with c GG /Λ. This is true also for the the relic line but only for m χ < m S , since for larger DM masses annihilation to mediators dominate and the line is effectively independent on c GG /Λ.
Unlike the photon fusion case discussed above, a thermal relic with pseudo-scalar mediator (right panel of Fig. 6) is less constrained in the gluon fusion regime. Limits from γ-ray lines searches are not applicable in this case, since the annihilation cross section in gluons (i.e. the one responsible for a continuum spectrum of photons) is up to 200 times bigger than the one in lines. Fermi limits from γ-ray continuum are of course still valid, but they exclude regions way above the thermal relic line. As for the scalar case, mono-jet searches do not put bounds in this DM mass range. In this case, the relic line is very smooth and the drop around m χ = m P is not visible.

DM Dominated Resonance
The other half of the parameter space corresponds to DM masses below m S /2, yielding the DM dominated resonance scenario discussed in Section 4.3. We use again the output of our LHC study to identify interesting classes of DM models.
We now examine the DM phenomenology in the photon fusion regime. Results are shown Figure 7: DM analysis in the DM dominated resonance scenario at the LHC discussed in Section 4.3 for scalar (left) and pseudo-scalar (right) mediators. In this figure we consider the photon fusion regime, with Wilson coefficients c BB in the range shown. We present results in the (m χ , c χ{S,P } c BB /Λ) plane, since rates at DM searches only depend on these two quantities. This is not exact for the thermal line due to resonance effects, and the green band around the thermal relic line quantifies how much this rescaling is violated. We shade away the region where the combination on the vertical axis is above 4π TeV −1 , and shade with gray where Γ S m S . The light gray region on the right of each figure corresponds to the SM dominated scenario already discussed. We show where thermal freeze-out can reproduce the observed DM density and shade regions excluded by DM searches. For Fermi limits, solid and dotted lines are for bounds coming from diffuse photons and lines, respectively. We identify the line reproducing the ATLAS preferred value for the total width Γ S 45 GeV. See text for further discussion.
in Fig. 7 for the case where only the coupling c BB is switched on. We present our results in a slightly different way here, putting the combination c χ{S,P } c BB /Λ on the vertical axis. Any rate for current and future experiments only depends on this combination, and therefore the same holds for the exclusion regions in the figures. However, it is less obvious that resonance effects on the relic density calculation [93] have a small impact, since the resonance is quite broad with a total width Γ S /m S c 2 χS /(8π), as follows from Eq. (37). We address this issue in each case and show that this effect is very small, hence not a concern for our final results.
We start the discussion from the left panel of Fig. 7, corresponding to a scalar mediator in the photon fusion regime. DM searches are powerless for this case 5 . The solid relic density line is explicitly obtained for Γ S 45 GeV. As shown in the top-left panels of Figs. 3 and 4, this corresponds to c BB /Λ 0.26 TeV, and this allows us to identify the iso-contour Γ S 45 GeV in 5 To give an idea, for 30 GeV DM and c BB /Λ 0.26 TeV −1 we find a DM-nucleon cross section of σ SI = 6.8 c 2 χS ×10 −53 cm 2 , well below the expected LZ sensitivity for c χS < 4π. As can be seen from Fig. 4, in principle c W W = c BB /Λ 0.26 TeV −1 is not excluded and in this case LZ could probe c χS values of O (10). Clearly, the chances of DD in the photon fusion regime are extremely slim. the (m χ , c χS c BB /Λ) plane of Fig. 7. A thermal relic consistent with Γ S 45 GeV then requires c χS 2.42. If resonance effects are negligible, the relic density line is a universal function of c χS c BB /Λ and we have explicitly checked that this rescaling invariance works perfectly for lower values of c χS . We expect it to break down for large enough c χS , and we estimate the error we could make with this rescaling by computing self-consistently the relic density line for a thermal relic with Γ S m S . As can be seen from Figs. 3 and 4, this corresponds to a larger coupling c BB /Λ 0.53 TeV, and a thermal relic would then require c χS 6.77. The net result on the relic density is a combination of two effects: a large overall coupling in the cross section and a broader width of the mediator. The green bands in the figure show that these combined effects are rather mild. Given the lack of constraints from DM searches, a scalar portal in the photon fusion regime leads again to a viable DM candidate. To summarize this discussion we present here the values of the parameters consistent with a thermal relic and a scalar width of 45 GeV:

�� ��
The LHC analysis allows for values of c W W c BB , but turning on this coupling does not greatly impact the obtained DM parameters m χ and c χS . The pseudo-scalar case is shown in the right panel of Fig. 7 with identical conventions. We again give the values of the parameters for a thermal relic and a 45 GeV width: However, in contrast to the scalar case, ID limits are now quite severe, and the γ-ray line bounds completely rule out a thermal relic even for a cored DM profile like the isothermal one. Finally, in Fig. 8 we consider the gluon fusion regime for a DM dominated resonance. In both panels, c BB /Λ and c GG /Λ are understood to be within the range written in the label, consistently with what found in our LHC study for the gluon fusion regime (see Figs. 3 and 4). We present our results in terms of the combination of couplings c χS c GG /Λ. DD cross sections only depend on this combination. The same is true for the relic density line, which we wish to draw again for Γ S 45 GeV. However, this choice does not identify the value of couplings to SM gauge bosons because the LHC analysis only fixes the product c GG c BB , as can be seen from the top-left plot of Figs. 3 and 4 or directly in Eq. (39). We therefore take the smallest value of the gluon coupling which is still inside the gluon-dominated regime, c GG /Λ 0.03 TeV −1 , which then yields c BB /Λ 0.01 TeV −1 . We see that LUX constraints are already very close for this gluon coupling and larger values (the LHC di-jet limit is c GG /Λ < 0.17 TeV −1 ) are in conflict with the bounds. Alternatively, keeping c GG /Λ 0.03 TeV −1 fixed, we see that Γ S /m S cannot be much larger than 6%. Lastly, mono-jet searches only put constraints on the Wilson coefficient c GG /Λ in this regime (see Eq. (40)). We show collider limits in this plane by choosing the value c χS = 2.91, giving a Γ S = 45 GeV width for a thermal relic and c GG /Λ 0.03 TeV −1 . Collider bounds are superseded by LUX at small masses, and they become relevant near the resonance. We conclude by giving explicit parameters for a benchmark point not excluded by LUX and LHC, consistent with a thermal relic and yielding the ATLAS preferred width: These values are in excellent agreement with those found in Ref. [11] and give a spin-independent DM-nucleon cross section σ SI 1.88 · 10 −45 cm 2 that can be probed in the near future by the LUX experiment (the current bound is 2.6 · 10 −45 cm 2 for m χ 310 GeV [79]). They could also yield a mono-jet signal at the LHC Run 2. As before, nonzero values of c W W c BB are not excluded but do not greatly impact the DM phenomenology.
The right panel shows the pseudo-scalar case. In this scenario, the ID limits from γ-ray lines depend on the specific value of c BB which is not fixed unlike the analysis in Fig. 6. We present our limits for c BB /Λ = 0.01 TeV −1 , which reproduces the ATLAS preferred width if one chooses the smallest allowed coupling to gluons in the gluon-fusion regime. The values of the parameters for a thermal relic and a 45 GeV width read: m χ 268 GeV , c χP 1.47 , c GG /Λ 0.03 TeV −1 , c BB /Λ 0.01 TeV −1 , c W W = 0 . (51) As one can see from Fig. 8, ID limits are very stringent and the γ-ray line bounds rule out a thermal relic. However, the LHC analysis only fixes the value of c GG c BB in this gluon fusion regime, unlike the photon-fusion case with parameters for DM given in Eq. (49). Thus we can choose a larger c GG and a smaller c BB , which of course makes the Fermi γ-ray line bounds less stringent. On the other hand, a larger value of c GG will also move up the 45 GeV width line in the (m χ , c χP c GG /Λ) plane by the same factor as the γ-ray line constraints, crossing the relic density line for smaller DM mass. We explicitly checked that a 45 GeV width can be obtained while evading the Fermi bounds for a DM mass of roughly 220 GeV and coupling to gluons c GG /Λ 0.05 TeV −1 . Here is where mono-jet bounds come into play. In analogy to what we have done for the scalar case, we present j+MET limits for c χP 1.47 as in the benchmark point of Eq. (51). We observe a different shape of the mono-jet bound as compared to the scalar case, consequence of the fact that the width of the resonance has a different dependence on the DM mass for scalar and pseudo-scalar. In particular, the j+MET limits are less stringent around the resonance for the pseudo-scalar case, because for fixed DM mass m χ 375 GeV and same couplings the decay width of the pseudo-scalar is typically larger. Nevertheless, this limits how much we can increase the coupling to gluons and therefore it makes a DM candidate with pseudo-scalar portal in the gluon fusion regime quite unlikely.

Outlook
In this paper we have studied the minimal EFT for the diphoton events recently observed at the LHC and DM. The field content is the same as the SM one with the addition of two gauge singlets, a (pseudo-)scalar and a Dirac fermion. We coupled the two singlets via a portal Yukawa interactions, and we also coupled the (pseudo-)scalar to SM gauge bosons via dimension 5 contact interactions. Due to observed decays to two photons, a coupling of the (pseudo-)scalar to electroweak gauge bosons is mandatory. On the contrary, the coupling to gluons is optional, as the new scalar can be produced through photon fusion in proton-proton collisions.
The LHC phenomenology turns out to be identical for scalar and pseudo-scalar, and we presented a study valid in both cases in Section 4. The knowledge of the resonance mass splits in two the possible values of the DM mass, according to whether invisible decays are kinematically accessible. We dubbed these two scenarios SM and DM dominated resonance, corresponding to DM masses that make invisible decays forbidden and allowed, respectively. Despite the six free EFT parameters (after fixing the resonance mass to 750 GeV), the parameter space region consistent with both the diphoton excess and bounds from LHC Run 1 are compactly summarized in Figs. 1-4. Remarkably, the Wilson coefficients are quite constrained and either gluon or photon fusion dominates the total production, unless we choose very specific ratios of the couplings. In the SM dominated scenarios we typically find a very small width for the new resonance. We have not attempted to construct explicit UV completions realizing the parameter space configuration identified by our analysis, consistently with the EFT spirit of this work. Explicit models in the gluon fusion regime have been studied in Refs. [24-31, 37, 94-102], and we think it would be very interesting to find some explicit realization of the photon fusion regime as well. Considering the large coefficients, the photon fusion scenario probably requires a strongly-coupled UV completion, see for instance Refs. [40,103]. Every sensible UV completion with a cutoff Λ few TeV should return Wilson coefficients at the LHC scale within the bounds identified in Figs. 1-4.
The results presented in Section 4 have a range of validity beyond DM models. Even in what we call the DM dominated case, our only assumption is the presence of some additional decay channels that does not have to be to neither stable nor cosmically abundant particles. But other than being interesting by itself, it significantly simplified our DM analysis in Section 5. We found it convenient again to split the DM discussion for the two different scenarios of SM and DM dominated resonance. Moreover, the two cases of scalar and pseudo-scalar mediator lead to drastically different DM phenomenology. Our findings can be compactly summarized by the following four classes of DM models: Scalar with Photon Fusion: DM searches cannot probe this parameter space region. Mono-jet bounds do not apply and ID rates are p-wave suppressed. The only hope would be direct searches, but the RG induced couplings given in Section 3.3 are below the LZ projected sensitivity. A thermal relic for DM masses above m S /2 but below m S can only be attained for DM couplings close to the perturbative limit, while for larger DM masses perturbative values of c χS are allowed. On the other hand, a thermal relic is totally viable for DM masses below the resonance, and the ATLAS preferred value Γ S /m S 6% can be achieved with invisible decays. Results are summarized on the left panels of Fig. 5

and 7.
Scalar with Gluon Fusion: DM annihilations are still p-wave suppressed, but LHC and DD experiments can put strong constraints. LUX bounds, evaluated through the RG prescription given in Section 3.2, are typically stronger than the ones from mono-jet. The only exception is for DM masses right below the resonant value m S /2, where LHC limits slightly overtakes the ones from DD. In this mass region a thermal relic is consistent with collider and direct searches, and it would give a signal in future experiments. On the other hand, LUX rules out most of the parameter space for masses between m S /2 and m S , while for masses above m S the parameter space is currently viable. This entire scenario, regardless of the specific value of the DM mass, will be deeply probed by LZ. Results are summarized on the left panels of Fig. 6 and 8.
Pseudo-scalar with Photon Fusion: DM annihilations mediated by a pseudo-scalar particles are s-wave processes. ID constraints are the only meaningful ones in this case, since DD rates are very suppressed and mono-jet limits do not apply. For m χ 500 GeV, Fermi searches for photon lines basically rule out a thermal relic. For larger masses the implications of H.E.S.S. limits are less obvious as they are quite sensitive to the density profile assumption. Results are summarized on the right panels of Fig. 5 and 7.
Pseudo-scalar with Gluon Fusion: Introducing a pseudo-scalar coupling to gluons has two main effects on DM phenomenology: making mono-jet searches meaningful and contaminating the line signals with consequent weakening of the ID constraints. Neither of these bounds quite gets to freeze-out line for DM masses above the resonance. The situation is rather different for DM masses smaller than m S /2, where the combination of limits from Fermi γ-ray line searches and LHC mono-jet searches is strong enough to rule out a thermal relic. Results are summarized on the right panels of Fig. 6 and 8.
If the diphoton excess turns out to be more than just a statistical fluctuation, we may have started a new era of discoveries in particle physics. Among other things, such as being part of a theoretical construct that solves the gauge hierarchy problem, this new particle could be the connector between the SM and the dark sector. Our general EFT analysis identified a broad class of DM models with a 750 GeV (pseudo-)scalar portal consistent with current experimental limits. Although the study of a specific DM theory goes beyond the purpose of this work, our results in Section 5 clearly pinpoints preferred models. The most appealing one is presumably the scalar mediator case in the gluon fusion regime, since it could soon give a signal in direct and collider searches. Contrarily, scalar portals in the photon fusion regime are unattainable by all DM experiments. In these cases, as well as pseudo-scalar cases in the gluon fusion regime and for large DM masses where ID limits are not as powerful, the most promising strategy to probe the models is to accumulate more evidence through other LHC channels such as Zγ searches. We believe it would be very interesting to further investigate the associated phenomenology of specific UV-complete DM models reproducing our EFT framework at lower energies.

Acknowledgments
We thank Lawrence Hall, Gordan Krnjaic, Marco Nardecchia, Duccio Pappadopulo, Diego Redigolo, Alessandro Strumia, Tim Tait for useful discussions. This work was supported by the U.S. Department of Energy grant number DE-SC0010107 (FD) and the Dutch Organization for Scientific Research (NWO) through a VENI grant (JdV). P.P. acknowledges support of the European Research Council project 267117 hosted by Université Pierre et Marie Curie-Paris 6, PI J. Silk.

A Decay Widths and Cross Sections
In this Appendix we give all the details of the results for decay widths and cross section for both LHC production and DM annihilation.

Interactions for Mass Eigenstates
In this first Appendix we express the interactions in Eq. (4) in terms of the mass eigenstates, and provide all the squared matrix elements for decays of the scalar S to any possible final state. The results contained here will be the building blocks to easily obtain decay widths and cross sections.
The SM charged EW bosons are obtained by a π/4 rotations among the SU (2) L gauge bosons The neutral gauge bosons are expressed by a weak mixing angle rotation The interactions in Eq. (4) can be rewritten as a function of the mass eigenstates with Wilson coefficients connected to the gauge invariant ones as follows The couplings for the pseudo-scalar case are analogous with just the replacement c XX →c XX , where X = {G, W, B} q a k 2 k 1 b i j Figure 9: Left: Feynman diagram for the S decay process to a generic two-body final state. In this Appendix we give the squared matrix elements for arbitrary external four momentum q of the scalar, and we define s = q 2 . The on-shell case corresponds to s = m 2 S . Right: Feynman diagrams for arbitrary annihilations 2 → 2.

Squared Matrix Elements
We use the above Lagrangian to compute the squared matrix elements for the decays process S → ij. The Feynman diagram is shown on the left of Fig. 9. We keep the external scalar state off-shell, and we define its invariant mass to be q 2 = s. In the on-shell limit, which we will take for example to compute decay widths, we will have s = m 2 S . We also sum over all the possible final polarizations. All the following calculation have been performed by hand and cross-checked with FeynCalc [104]. Denoting by k 1 and k 2 the four momenta of the final state particles, we find the following squared matrix elements for decay processes to SM final states Likewise, we can evaluate the squared matrix element for decay to DM pairs We switch to the case of a pseudo-scalar, with matrix elements for decays to SM states Finally, the squared matrix element for decays to DM results in

Decays Rates
With the squared matrix elements in hands, it is straightforward to compute the partial decay width for a generic channel. We have the general expression for scalar decays where the statistical factor s ij accounts for identical particles in the final state. We find The expression for pseudo-scalar decays can be obtained identically.

LHC Cross Sections
The total cross section for the process pp → ij is obtained from the factorization theorem where f a/p and f b/p are the a and b parton distribution function (pdf) inside the proton. We introduce a new final state variable x, defined as the invariant mass of the ij-pair in units of m S . By using energy-momentum conservation we can rewrite this variable as follows The total cross section can written in a compact form where we define the flux function at fixed CM energy for the pp collision The partonic cross section for a 2 → 2 collisions takes the general form where as usual s ij accounts for identical particles in the final state and we define the function g(x, y) = 1 − 2(x + y) + (x − y) 2 .
The above equations are general. We specialize now to the case of a s-channel resonance shown in Fig. 9, where the partonic matrix elements always take the form The expression for a s-channel pseudo-scalar resonance P is identical. Here, the factor of 1/4 average over the initial polarizations, since every possible initial state has always 2 polarizations. We also account for a possible average over the color number N ab , as we can have gluons in the initial state. Furthermore, we specialize to the case of only gluons and photons in the initial state, and we write the partonic cross section The partial decay widths in the above equation have to be computed as we would for a scalar particle S of mass √ s.
The invariant mass of the diphoton pairs observed at the LHC is never far from m S , therefore we can further simplify our expression by employing the narrow width approximation The dx integral in Eq. (80) is straighforward where we define

DM Annihilation Cross Sections I: SM Final States
We collect the total cross sections for DM annihilation to SM final states. Considering the process χχ → ij, the cross section formally reads Here, √ s is the energy in the CM frame of the collision and the statistical factor s ij = 1/2 for identical particles in the final state.
The squared matrix element for the collision mediated by a scalar exchanged in the s-channel can be expressed as follows (see Fig. 9) where the factor 1/4 averages over the initial DM polarizations. Plugging the squared matrix element for S → ij as given in Eq. (64), we find the general expression for the DM annihilation cross section The result for each channel ij can be found by plugging the squared matrix elements given in this Appendix. Likewise, the expression for processes mediated by a pseudo-scalar results in For the last two equations, we see that annihilations mediated by the scalar and the pseudoscalar are p-and s-wave processes, respectively.

DM Annihilation Cross Sections II: Mediators Final States
DM annihilations to mediator particles become kinematically accessible for m χ 750 GeV. These processes are mediated by a virtual DM particle exchanged in both t-and u-channels. We computed the full cross section as a function of the CM energy √ s and used them for the relic density calculation. The general expressions are quite involved. In this Appendix we only report the non-relativistic limits for annihilation to scalar and pseudo-scalars where The processes are p-wave suppressed in both cases. These approximate results are quite accurate since we are away from the resonant value m χ 375 GeV, but nevertheless we use the full expressions for our numerical analysis.

B RGE: Equations and Solutions
In this Appendix we give the details of our RG analysis. As done in Section 3, we divide the discussion into two cases according to whether we have a coupling to gluons at the cutoff scale.

RGE with coupling to gluons at the cutoff
For non zero values of c GG (Λ) and in the renormalization scale range m S < µ < Λ we limit ourselves to the following effective Lagrangian with H the SM Higgs doublet and y q the quark Yukawa couplings. For this RGE analysis we find it convenient to employ a different normalization for the coupling to gluons c GG (µ)α s (µ) = c GG (µ) , as it yields a simple anomalous dimension matrix. Likewise, factorizing out the Yukawa couplings ensures that the evolution of c yq is the same for every quark. We define the two dimensional array of Wilson coefficients c(µ) T = (c yq (µ) c GG (µ)) and write the RG equation as d c(µ) d ln µ = γ(µ) c(µ) , in terms of the anomalous dimension matrix γ. As anticipated, the normalization chosen in Eq. (97) ensures the simple form where γ Gq describes the mixing from c GG into c yq . We work in a mass independent renormalization scheme (such as MS) where the anomalous dimension matrix depends on the renormalization scale µ only through SM couplings. Throughout this work we use leading-order (LO) QCD evolution such that [105] γ Gq (µ) = α 2 s (µ) 4π where γ 0 = 32. We always take c yq (Λ) = 0 as our boundary condition. The solution of the RG system in the energy range m S < µ < Λ allows us to obtain the couplings at the scale µ = m S . We use the output of this operation in Section 3.1 to justify how such a running has a negligible impact on LHC phenomenology. For this reason, we neglect this running and start our actual RG analysis at the scale µ = m S , with the coupling c yq still set to vanish. At this scale, we integrate out the scalar resonance and obtain the effective Lagrangian with boundary conditions C yq (m S ) 0 , With this choice of boundary conditions we start our RG evolution at µ = m S and connect the couplings in Eqs (102) and (103) with the ones at the nuclear scale. The evolution down to the electroweak scale, where we integrate out the top quark, goes along almost the same exact lines. The main difference is that a threshold correction to C GG is induced after the top quark is integrated out [106]. At slightly lower energies, we break the electroweak SU (2) L × U Y (1) → U (1) via the Higgs vacuum expectation value. For simplicity, we perform these steps at the same scale m t and we match to the effective Lagrangian where, at a scale right below µ = m t , C q (m t ) = C yq (m t ) = γ 0 2β 0 (α s (m S ) − α s (m t )) C GG (m S ) , C GG (m t ) = C GG (m S ) − 1 12π C yt (m t ) .
Here, β 0 = 11 − 2n f /3 and n f is the number of active flavors (n f = 6 for µ > m t ).
The evolution to µ N is now straightforward, with the main difference being that the number of active quark flavors is reduced by one after each quark threshold. We give here the numerical result for the evolution from m S to µ N which is independent of the cut-off scale Λ. For q = {u, d, s}, we obtain C q (µ N ) = − 0.54 C GG (m S ) , (107) C GG (µ N ) = 1.02 C GG (m S ) .
We note that the 2% correction in Eq. (108) is actually a two-loop effect as it arises from C GG mixing into C q and a subsequent threshold correction. Its small size indicates that the LO analysis is sufficient for our purposes. Using the values α s (m S ) = 0.092 and α s (µ N ) = 0.362, we can derive the low-energy couplings in the language of the basis of Eq. (7), as given in Eqs. (16) and (17) of Section 3.

RGE without coupling to gluons at the cutoff
In the second scenario we assume that S does not couple to gluons at the cutoff scale Λ. That is, the effective Lagrangian at the scale µ = Λ has only the interactions The RGE to lower energies requires the inclusion of two additional operators [20] ∆L EFT = S Λ q=u,d,s,c,b,t c yq y q (q L Hq R + h.c.) + c H (H † H) 2 .
Also in this case we take c yq (Λ) = c H (Λ) = 0 as our boundary conditions. However, they are radiatively induced at lower scales and must be kept in order to have a consistent RG analysis. Analogous to the previous scenario, we define c(µ) T = (c yq (µ) c H (µ) c BB (µ) c W W (µ)) and write the RG equation as d c(µ) d ln µ = γ(µ) c(µ) , in terms of the anomalous dimension matrix γ(µ). As we are now considering electroweak corrections, we only consider the mixing of c BB and c W W into c yq and c H . We neglect the evolution of g, g , c BB and c W W themselves as this would correspond to α 2 em corrections to direct detection cross sections. That is, we approximate The anomalous dimensions can be straightforwardly calculated and we find 6 γ BH = − 6 (4π) 2 g 4 + g 2 g 2 , γ W H = − 6 (4π) 2 3g 4 + g 2 g 2 .
where Q q denotes the quark charge in units of |e| (e.g. Q t = +2/3). At the scale m S we integrate out the scalar resonance and work with the effective Lagrangian L mt<µ<m S EFT = C W Wχ χ W I µν W I µν + C BBχ χ B µν B µν + , q=u,d,s,c,b,t C yq y qχ χ (q L Hq R + h.c.) + C H (H † H) 2χ χ .
As done for the QCD running, we set our boundary conditions at the scale µ = m S C yg (m S ) 0 , and do not account for the negligible running from Λ to m S . We do evolve the couplings from m S to the scale of electroweak symmetry breaking. As before, we integrate out the heavy SM degrees of freedom at the common scale m t and then match to the effective Lagrangian C q m qχ χ qq + C GGχ χ G A µν G A µν + C F Fχ χ F µν F µν .
We find the following Wilson coefficients where m + t denotes a scale slightly above the top quark mass and v = 246 GeV. Note that here we neglected an O(α 2 em ) threshold correction to C F F from integrating out the top quark.
We can evolve this set of operators to lower energies. As C GG (m t ) is only induced at the two-loop level, we neglect additional mixing from C GG into C q . The only mixing we then need to consider is the mixing between C F F and C q described by the anomalous dimension [14] γ F q = α em 4π (24Q 2 q ) .
The factor of 2 in Eq. (134) is because we deal with a Dirac fermion and we add the contribution of the antiparticles. Finally, we compute the DM contribution to the Ω parameter where for the critical density we have [107] ρ cr /h 2 = 1.05375 × 10 −5 GeV cm −3 .