The Dark Penguin Shines Light at Colliders

Collider experiments are one of the most promising ways to constrain Dark Matter (DM) interactions. For several types of DM-Standard Model couplings, a meaningful interpretation of the results requires to go beyond effective field theory, considering simplified models with light mediators. This is especially important in the case of loop-mediated interactions. In this paper we perform the first simplified model study of the magnetic dipole interacting DM, by including the one-loop momentum-dependent form factors that mediate the coupling -- given by the Dark Penguin -- in collider processes. We compute bounds from the monojet, monophoton, and diphoton searches at the $8$ and $14$ TeV LHC, and compare the results to those of direct and indirect detection experiments. Future searches at the $100$ TeV hadron collider and at the ILC are also addressed. We find that the optimal search strategy requires loose cuts on the missing transverse energy, to capture the enhancement of the form factors near the threshold for on-shell production of the mediators. We consider both minimal models and models where an additional state beyond the DM is accessible. In the latter case, under the assumption of anarchic flavor structure in the dark sector, the LHC monophoton and diphoton searches will be able to set much stronger bounds than in the minimal scenario. A determination of the mass of the heavier dark fermion might be feasible using the $M_{T2}$ variable. In addition, if the Dark Penguin flavor structure is almost aligned with that of the DM mass, a displaced signal from the decay of the heavier dark fermion into the DM and photon can be observed. This allows us to set constraints on the mixings and couplings of the model from an existing search for non-pointing photons.


Introduction
The existence of Dark Matter (DM) is firmly established by a large number of astrophysical and cosmological observations. Despite the fact that it contributes a large component of the energy density of the universe, however, its precise properties remain almost completely mysterious. The common belief is that most of the DM is in the form of a stable particle, which is neutral or charged very weakly under the electric force, but interacts at least gravitationally with baryons. If such a particle carries non-gravitational interactions with the Standard Model (SM), the production of DM particles at high energy collider experiments offers one of the most promising opportunities to identify the nature of the DM interactions. Thanks to the capability to produce DM particles in a wide mass range and to the obvious independence from astrophysical uncertainties, collider searches provide complementary results to the direct and indirect detection experiments [1][2][3][4][5][6].
When setting collider constraints on DM interactions, the most straightforward way of parameterizing the DM-SM coupling is through an effective field theory (EFT) description. The non-observation of events with significant missing energy in excess of the SM background is then translated into upper bounds on the coefficients of the effective operators that couple the DM to the SM fields [7][8][9][10][11][12][13][14][15][16][17][18][19][20][21]. This simple method is independent of the details of the ultraviolet (UV) completion of the model. However, it is based on the assumption that the EFT gives a valid description of the collider process, i.e. that the mediators can be integrated out at the energy scale of the collision. Unfortunately, the sensitivity of the current LHC searches does not correspond to heavy enough mediators for many of the EFT couplings [22,23], and unitarity usually sets stronger bounds than the collider search [24,25]. Thus in many instances, to extract meaningful information from collider searches it is necessary to consider perturbative models with light mediators. This 'simplified model' approach has been applied to several DM-SM couplings whose UV completions feature the tree-level exchange of mediators [26][27][28][29][30][31][32][33][34][35][36][37][38][39].
The inadequacy of the EFT approach is manifestly even more dramatic in the case of loopmediated couplings. Among the loop-induced DM-SM interactions, an especially important status is held by the dipole and Rayleigh operators [40,41], which play important roles in DM model building. Collider constraints on these couplings within the EFT approach have been explored in several previous studies [42][43][44][45][46], but the resulting bounds translate into very weak constraints on the mediator masses, thus calling for a simplified model description.
In this work we perform the first systematic collider study of a loop-induced DM-SM coupling, by considering a perturbative UV completion with light mediators and including the momentum-dependent form factors in the description of collider processes. 1 In our study we focus on a generic UV completion of the magnetic and electric dipole operators, the Dark Penguin. We emphasize, however, that the method introduced in this work can be applied to any other loop-induced DM-SM coupling.
If the DM is a singlet under the SM gauge symmetry, the magnetic and electric dipole operators µ Mχ σ µν χF µν + µ Eχ σ µν χF µν , generated by the dark penguin diagrams, give the lowest order interactions of the DM with the SM gauge fields. As a consequence, these operators play important roles in the possible explanation of various γ-ray excesses observed at indirect detection experiments [49][50][51][52][53][54][55]. Moreover, being chirality-flipping, these operators 1 A first brief discussion of the loop form factors for dipole-interacting DM production at colliders was presented in Ref. [47]. See also Ref. [48], where the loop contribution of the dark sector to dilepton production was considered, although in a different model.
can possess a non-trivial flavor structure when more than one species of dark fermions χ i is present. In this case the photon dipole couplings in the mass basis can connect the light and excited DM state, and give interesting inelastic scattering signals at direct detection experiments [56][57][58][59][60][61]. The typical momentum exchange in the direct and indirect detection experiments is very small, therefore it is reasonable to use the EFT description even when the mediators are light with respect to collider energies.
In order to allow a comparison to the results of direct and indirect detection experiments, the collider searches need to provide bounds on both the DM-mediator coupling and the mediator mass. The search for mediator decays in various UV completions of the dipole and Raleigh operators has been studied in [62]. However, since the decay of the mediators only depends on their branching ratios, the information about the DM-mediator coupling is lost in these searches. Thus the study of the dark penguin process would be crucial even if the mediators were directly discovered first.
The dark penguin serves as a good example of a simplified model that can be constrained by different collider searches: depending on the 'dark flavor' structure -the flavor structure of dark fermions, -it can provide signals in the monojet, monophoton, diphoton, and even non-pointing photon searches. Besides allowing us to set meaningful constraints, the inclusion of light mediators helps us to identify the optimal cuts to be used in the collider searches, which are different from those commonly employed in the study of the EFT couplings. When the mediator mass is much smaller than the typical momentum exchange in the dipole, the missing transverse energy (MET) distribution of the DM signal is much softer than the one obtained assuming EFT couplings. Therefore the intuition of setting harder cuts to increase the signal excess no longer applies: on the contrary, the best strategy is to keep the cuts as low as possible, to include the enhancement of the form factors corresponding to the mediators being produced on-shell. It follows that in the search for the dark penguin with light mediators, lowering the background is more important than increasing the collider energy. Looking ahead towards future experiments, this implies that the International Linear Collider (ILC) would have the capability to set much stronger constraints than a very high energy (100 TeV) hadron collider, as will be shown in our analysis.
The organization of the paper is as follows. In Sec. 2 we present the simplified model used in the analysis, and address the possibility of having displaced photon signals from an aligned flavor structure. We begin Sec. 3 by explaining our method for including the loopmediated dipole couplings in the simulation of collider processes. Then we discuss in detail the missing energy searches used to set constraints on the simplified model, paying particular attention to the estimation of the systematic uncertainties in the projection to the 14 TeV LHC. In Sec. 4 we present the constraints from LEP, the 8 and 14 TeV LHC, and the future 100 TeV collider and ILC on the dark penguin parameter space, assuming a single flavor of dark fermions. Furthermore, we point out that the sensitivity to light mediators can be improved by choosing MET cuts weaker than those used in the search for effective couplings.
To give an example of the case with more than one dark flavors, we analyze the model with two dark fermions in Sec. 5, setting bounds from LEP and the monophoton and diphoton searches at the LHC. For the diphoton channel, the possible application of the M T 2 variable to determine the mass of the heavier dark fermion is discussed. In Sec. 6 we consider the dark penguin with a flavor structure almost aligned to the dark fermion masses, in which case displaced photon signals can be observed at the LHC. We use the 8 TeV ATLAS search for non-pointing photons to compute the bound on the dark mixing angle and coupling. In Sec. 7 the collider constraints on the dark penguin are compared to the current results from direct and indirect DM detection experiments, by showing the reach of the different searches on the magnetic dipole moment and the annihilation cross section of χχ → γγ, respectively. We conclude by summarizing our result in Sec. 8. Finally, App. A contains general formulas for the dark penguin form factors, whereas App. B collects the basic statistics we used for setting limits.
It is worth pointing out that Sec. 3 is somewhat technical, therefore the reader mainly interested in the results of our work might prefer, after having become familiar with the dark penguin in Sec. 2, to move directly to Sec. 4.

A simplified model of dark penguin
Here we describe a simple UV completion of the magnetic and electric dipole operators, which will be employed throughout the paper. The model is similar to the one discussed in [49,62] and contains dark Dirac fermions χ i with flavor index i = 1, ..., N χ , a fermion mediator ψ, and a scalar mediator φ. Both mediators carry some hypercharge Y , and the χ i are SM singlets. A specific assignment of dark charge which stabilizes the DM particle may affect the decay of mediators, but not the dark penguin process we are interested in. The Lagrangian in the mass basis is written as Several simplifications will be made on the Lagrangian: we assume the fermion and scalar mediators to have the same mass, M f = M s = M , 2 and also assume the Yukawa-type couplings λ i = λ L i = λ R i to be real. This in particular implies that no electric dipole moment operator is generated from Eq. (2.1). Depending on the details of a specific model, one of the mediators can decay into SM particles, while the mediator that carries the dark charge stabilizing the DM can decay into the DM and SM particles. A general analysis of the phenomenology of the mediators can be found in [62]. Their results show that the discovery prospects at the LHC depend strongly on the charge assignment. The most challenging scenario corresponds to SU (2) L -singlet mediators with Y = −1, and ψ mixed with the right-handed tau lepton. Even at 14 TeV with 300 fb −1 , the LHC will have no sensitivity to this model [62]. However, different hypercharge assignments can change dramatically the decays of the mediators, leading to significantly better prospects. In our analysis of the dark penguin we wish to be independent from the details of the model building, therefore we Figure 1. Diagrams of the dark penguin. The corresponding amplitude is given in Eq. (2.2).
treat Y and M as free parameters. Precision electroweak measurements do not set relevant constraints on mediators carrying only hypercharge. Even if the mediators are doublets under SU (2) L , as considered for example in Ref. [49], no contribution to the S and T parameters arises at one loop.
The model in Eq. (2.1) generates a magnetic dipole operator through the dark penguin process in Fig. 1. The amplitude for B µ → χ iχj is written into a gauge invariant form In this paper, we use various collider searches to set an upper bound on 3 λ √ Y N .
where the factor λ i λ j encodes the dark flavor structure. The factor N is the multiplicity of the mediators (N > 1 can arise, for example, if the mediators transform non-trivially under SU (2) L ). The form factors F q,σ as functions of the masses (m χ i , m χ j , M f , M s ) and momentum q 2 (defined in Fig. 1) are given in App. A. It is important to note that the model discussed here can also generate other loops, such as those for the Rayleigh operator,χ i χ j B µν B µν . However, the contributions of this operator to the collider processes considered in this paper carry either extra gauge couplings or phase space suppressions, and thus only give sub-leading effects, except in a few cases on which we will comment in what follows. We leave the detailed study of the Rayleigh operators for future work.
In presence of a non-trivial flavor structure of the dark fermions χ i , one important difference between collider searches and (in)direct detection experiments is that the collider processes can generically involve more than one dark flavor. Indeed, as we will show, the current and upcoming missing energy searches at colliders are strikingly more sensitive to the case where more than one dark flavors are within kinematic reach. In this case, some assumptions on the flavor structure are necessary in order to compare results between the different experiments. In this paper, we assume the λ i couplings have a totally anarchic structure, with no unnatural hierarchies between different flavors. The assumption of anarchic structure permits a direct comparison of the collider bounds among themselves and to the direct and indirect detection experiments, which are typically sensitive only to the couplings of the lightest dark fermion (the DM particle). To be more precise, in the following we assume For models with different flavor structures, one can rescale our bounds from the different searches to set proper constraints.
According to naive dimensional analysis (NDA), the perturbative bound is λ √ N ∼ < 4π. When showing our results, we allow λ √ Y N to be as large as 8π, which corresponds to the NDA perturbative bound for Y = 4 (corresponding to g Y 1). In a different UV completion of the dipole interaction, such as for example a model where the scalar mediator is replaced with a gauge boson, the quantitative result for the dark penguin would be modified. However, the fact that gauge invariance constrains the structure of the amplitude forbids a qualitative alteration of our analysis.

Displaced signals and the aligned flavor structure
Another interesting scenario to explore is the case where the flavor structures of the DM mass matrix and the dark penguin are nearly aligned. This means that the fermion mass matrix m ijχi χ j is almost proportional to the dark penguin, whose flavor structure is determined by where the term proportional to the identity gives the lightest dark fermion a mass ∼m. When rotating the dark fermion fields into the mass basis, the dark penguin between different mass eigenstates carries an extra O( ) suppression that can result into a displaced decay of a heavy dark fermion χ h to the light DM χ l and a photon, with a width (assuming M m χ h m χ l ) Contrarily to the case where the decay is displaced because of phase space suppression [63], here the photon is hard, therefore the process is accessible in LHC searches. The aligned flavor structure can be generated if the dark sector has a single flavor breaking spurion, and both the DM masses and dipole interactions are generated by the same loop-level mediation. Since the gauge coupling is flavor blind, the flavor structures of the two operators are identical implying that no heavy-light dipole coupling is present in the mass basis for dark fermions. However, if there exists a small chiral symmetry breaking in the infrared, or an extra flavor symmetry breaking from an even higher order mediation is present, the misalignment between the newly generated mass and the dark penguin leads to a small heavy-light dipole coupling and thus a suppressed decay rate as in Eq. (2.5). With the assumption of anarchic couplings λ i , the search of non-pointing photons+MET can set interesting bounds on the (λ, ) plane. The result is discussed in Sec. 6.

DM searches at colliders
In this section we describe the collider analyses used in this paper to set bounds on the DM parameter space. We begin by outlining our procedure for taking into account the full loop dark penguin form factors in the collider simulations. Then we move on to describe in detail the relevant searches, including the monojet, monophoton and diphoton final states at the 8 and 14 TeV LHC. Then we turn to the 8 TeV search for non-pointing photons. Finally, the projections to a 100 TeV collider and to the ILC are studied for some of the searches. In our analysis of 8 TeV searches, we reproduce the shape of each SM background using our MonteCarlo (MC) simulations, and compute the additional global rescaling factor needed to obtain exact agreement with the distributions reported in the experimental papers. We then apply these rescaling factors in the 14 TeV projections. In the 14 TeV projection of monojet and monophoton searches, we estimate the improvement of systematic uncertainties by separating them into two parts. For uncertainties that relate to the normalization of the SM background, we follow the data-driven analysis at 8 TeV by simulating the control region sample at 14 TeV. For other uncertainties, including those on the PDFs and the acceptance of the detector, we assume an improvement proportional to the square-root of the luminosity. We do not attempt to simulate the QCD background at 14 TeV in this work. In the 8 TeV search, the pure QCD processes contribute less than ∼ 1% (10%) of the monojet (monophoton) background, and we make the reasonable assumption that the background can be kept subdominant at 14 TeV by setting a harder MET (photon p T ) cut. On the other hand, the QCD background does play an important role in the 8 TeV diphoton+MET search. We include additional jet and lepton vetoes to suppress it at 14 TeV.
We generated both the signals and backgrounds using MadGraph5 [64] and showered the parton level events using Pythia 6 [65]. We used PGS 4 for the detector simulation and cross-checked the results using Delphes 3 [66]. To compute the QCD K-factors for some SM backgrounds, we used MCFM [67] and VBFNLO [68].

Including loop form factors in a collider process
To properly describe loop-mediated processes at colliders, we simulate the DM production using EFT operators and reweight the events by employing the expressions of the form factors given in App. A. As described in Eq. (2.2), the amplitude of the dark penguin contains three distinct Lorentz structures The coefficient of the last operator ∼ q µ vanishes in the N χ = 1 case, since m χ i = m χ j . For N χ > 1 the contribution of this operator to the amplitude for pair production of DM particles ff → χ iχj is proportional tov f ( / p f + / pf )u f and is thus strongly suppressed by the light SM fermion masses, after using the equations of motion. Therefore we neglect the operator ∼ q µ altogether and only consider the σ µν q ν and q 2 γ µ operators.
For a given process, we first generate the events in MadGraph5 using a linear combination of the two relevant effective operators, then we reweight the events using the ratio of the partonic matrix element squared computed retaining the form factors to the value obtained using the EFT.
To give a concrete example, we consider the N χ = 2 monophoton analysis. We first simulate the signal process in Fig. 13 After performing the detector simulation, we apply the cuts for the monophoton search and obtain a list of signal events. For each event, we obtain the value of q 2 from the corresponding parton-level four momenta and use it to compute the form factors F σ,q (q 2 ). We then reweight the event using the expression of the matrix element squared, schematically where M Fσ,q (q 2 ) are the amplitudes corresponding to the first and second term in Eq. (3.1), respectively, whereas M O σ,V are the amplitudes corresponding to the effective operators. The procedure described above fully accounts for the interference between the two relevant Lorentz structures, and was applied in the study of the N χ = 2 scenario, where both monophoton and diphoton signals are given by 2 → 2 scatterings followed by the decay of χ h , as shown in Figs. 13 and 16. On the other hand, for N χ = 1 the main constraint comes from the monojet process in Fig. 10, which is genuinely a 2 → 3 scattering. For the sake of simplicity, in this case we neglect the interference term in Eq. (3.3), which carries an additional ∼ m χ / q 2 suppression due to the different chirality structure between the two operators. In the N χ = 2 case, the inclusion of the interference term gives at 14 TeV an increase of the cross section of less than 40%, which translates into a correction of less than 10% to the constraint on λ √ Y N . Such a deviation is acceptable when compared to the possible uncertainties in our projection. Thus we simulate the monojet events using the incoherent sum of the two operators in Eq. (3.2) to obtain the dark penguin result.
To demonstrate the effect of the form factors, we plot in Fig. 2 the q 2 distribution of the monojet events in both the dark penguin and EFT cases. Notice that due to the kinematics of the process, the MET is given by q 2 multiplied by some angular factors. The distribution of the dark penguin is remarkably different from the one of the EFT: the dark penguin exhibits an enhancement around q 2 2M , corresponding to the threshold for the production of on-shell mediators [47], and it is softer than the EFT at high energy. As a consequence, contrarily to the strategy used for effective DM couplings, where a strong cut on the MET is generically preferred, in the search for the dark penguin with light mediators a softer cut is favored. We will return to this point when discussing the monojet result, see Fig. 12. The preference for soft MET cuts applies also to the monophoton and diphoton searches.

LHC monojet
For the monojet channel, we follow and extend the CMS analysis in [70]. The event selection requires one jet with p T (j 1 ) > 110 GeV and |η(j 1 )| < 2.4. A second jet with p T (j 2 ) > 30 GeV and |η(j 2 )| < 4.5 is allowed, as long as the ∆φ(j 1 , j 2 ) < 2.5. Events containing a third jet satisfying p T (j 3 ) > 30 GeV and |η(j 3 )| < 4.5 are vetoed. To reduce the backgrounds from Z and W production, events containing electrons with p T > 10 GeV and |η| < 2.5, muons with p T > 10 GeV and |η| < 2.1, or hadronic taus with p T > 20 GeV and |η| < 2.3 are also vetoed. The counting experiments in [70] are performed in 7 signal regions, with MET cut from 250 to 550 GeV with a step of 50 GeV. In the left panel of Fig. 3 we show the comparison of our simulated backgrounds with those reported by CMS. For each of the three main backgrounds Z+jets, W +jets and tt, we use MCFM to obtain the QCD K-factor and apply an additional rescaling factor to match the CMS result. The rescaling factors are 0.84, 0.95 and 0.60, respectively. The rescaling factor for Z+jets is also applied to the DM signal events.
The CMS search uses a data driven analysis based on a µ+jets control sample to determine the normalization of the jZ and jW backgrounds. The corresponding events in the control   sample are Z(µ + µ − )+jets and W (µν)+jets, respectively. The cuts for the control sample are the same as for the monojet search, except the lepton vetoes are not applied. The Z(µ + µ − )+jets sample requires two muons with p T > 20 GeV and |η| < 2.1, with at least one of the muons passing isolation requirements, and the invariant mass of the muon pair between 60 and 120 GeV. Similarly, the W (µν)+jets sample requires one isolated muon with p T > 20 GeV and |η| < 2.1, and the transverse mass of the muon plus neutrino system in the range 50 GeV < M T < 100 GeV. Fig. 4 shows the comparison between our MC simulation of the control regions and the CMS results. For Z(µ + µ − )+jets, the MET is defined as the sum of the muon transverse momenta, while for W (µν)+jets the MET is given by the neutrino transverse momentum. The good agreement with the CMS results allows us to simulate the data driven analysis in the 14 TeV study.
For the 14 TeV projection we follow the same cuts in the CMS 8 TeV analysis, apart from  Table 1. Summary of the contributions (in %) to the uncertainty used for the 14 TeV monojet study, following the analysis in [70]. The dominant uncertainty comes from the limited number of control sample events.L ≡ L/(20 fb −1 ), and we assume the other sources of uncertainty will be improved with the increase of data.

CR statistics Other systematics
varying the MET cut from 550 to 2250 GeV with a step of 100 GeV. The dominant systematic uncertainty for each MET cut choice is obtained from the control sample simulation. For the remaining uncertainties, we take the values quoted in the CMS analysis for a MET cut of 550 GeV, and assume they will decrease with the square root of luminosity, see Table 1.
Since at 8 TeV these additional uncertainties do not depend strongly on the MET cut, it is a reasonable guess to use their value also at 14 TeV. The right panel of Fig. 3 shows the projected monojet backgrounds studied in this work, together with the dark penguin signal. The total systematic uncertainty is shown in Fig. 5 as a function of luminosity and MET cut. To provide some figure of merit, with 3 ab −1 of data and MET cut of 550 GeV the background is ∼ 10 6 events, which requires ∼ 10 4 DM signal events for a few-σ excess, given an uncertainty of order 1%. The corresponding λ √ Y N is ∼ 20 when the mediator mass is M 500 GeV.

LHC monophoton
For the monophoton channel, we follow the 8 TeV CMS search in [71]. The event selection requires one photon with p γ T > 145 GeV and |η γ | < 1.4442, and in addition / E T > 140 GeV with ∆φ( / E T , γ) > 2. Events containing electrons or muons with p T > 10 GeV, |η| < 2.5, and ∆R γ > 0.5 are vetoed, as well as events containing more than one reconstructed jet with p T > 30 GeV, |η| < 2.5 and ∆R j γ > 0.5. The comparison of our simulation of the photon p T distribution to the one obtained in the CMS analysis is given in Fig. 6. For each background we apply a p T -dependent K-factor obtained from the NLO MCFM calculation [67]. 5 In addition, we apply an overall rescaling in order to match the normalization provided by CMS: the rescaling factor is equal to 1.1 for Z(νν)γ, 1.2 for W (eν) and 1.8 for W ( ν)γ. 6 We do not attempt to simulate the backgrounds given by the misidentification of leptons (W (µν), Z(ll)γ), jets (γ j), γγ, beam halo, and QCD, but these backgrounds are subdominant as is shown in the left plot of Fig. 6. Only the three dominant backgrounds Z(νν)γ, W (eν), W ( ν)γ are included in the 14 TeV analysis.
In the CMS analysis [71], the main source of systematic uncertainty is given by the higherorder QCD corrections to the Z(νν)γ and W ( ν)γ backgrounds. Although the collaboration performs a data-driven analysis for each of the two backgrounds and compares the results to those obtained from MCFM, the limited number of control sample events at 8 TeV does not 5 The K-factor is very large for W ( ν)γ, K 4.5 in the region considered p γ T > 145 GeV and growing with increasing p γ T . This is a consequence of the presence of a 'radiation zero' [72,73] that strongly suppresses the LO amplitude. The NLO QCD corrections (in particular those associated with the qg partonic channel, which is not suppressed by the radiation zero) therefore contribute a large fraction of the total cross section, leading to the big K-factor [74]. 6 We do not apply to the DM signal the rescaling factor 1.1 obtained from the Z(νν)γ background, since it amounts to a negligible correction.  help to reduce the uncertainty. Therefore, in our 8 TeV analysis we made directly use of the uncertainties quoted by CMS. However, following the same strategy of the monojet case, we are going to assume that the data-driven analysis will play an important role in the future.
With the increase of luminosity, the systematic uncertainty will be based on the statistics of the control sample. In this case, using too stringent kinematic cuts in the search can increase the size of the systematic uncertainty, and a proper choice of the cut is required to optimize the signal significance.
To perform a data-driven analysis at 14 TeV, we compare our 8 TeV simulation of the two control region samples to the CMS results, and then apply the same analysis for the 14 TeV case. The Z(νν)γ control region is defined by the inverted angular cut ∆φ( / E T , γ) < 2.9, while the W ( ν)γ control region is defined by the inverted lepton veto. In both cases we obtain a result in good agreement with CMS, see Fig. 7. 7 This allows us to simulate the control region sample at 14 TeV and derive the relative systematic uncertainty shown in Fig. 5. As given in Table 2, the total uncertainty also contains various systematic uncertainties in the control region analysis for Zγ and W γ, as well as the uncertainty on the probability for an electron to be misidentified as photon for the W (eν) background. We take the values of these uncertainties from the 8 TeV analysis [71], and assume they will be improved with the increased luminosity as L −1/2 . This is based on the fact that various data-driven analyses have been used to determine these uncertainties at 8 TeV.
For the 14 TeV monophoton projection, we follow the same event selection used in the 8 TeV CMS search, except for tighter cuts on the photon p T and MET: both are required to be larger than 300 GeV, to suppress the QCD and beam halo backgrounds, which are not included in the 14 TeV estimate. As discussed in Sec. 3.1, choosing a tighter cut does not improve the sensitivity to the dark penguin. 7 Notice that we need to apply large rescaling factors for the W ( ν)γ and W (eν) backgrounds in the Zγ control region, 4.2 and 3.0, respectively. However, these backgrounds are subdominant to Z(νν)γ.
-13 -Background Relative systematic uncertainty Table 2. Summary of the contributions to the uncertainty used for the 14 TeV monophoton search. L = L/(20 fb −1 ), and the number of control sample events depends on the photon p T cut. The relative uncertainties are taken from the 8 TeV CMS analysis [71].  Table 3. Cuts, number of predicted events and number of observed events for the signal regions used in the 8 TeV ATLAS diphoton+MET analysis and in our 14 TeV projection. For the latter a luminosity of 300 fb −1 was assumed.

LHC diphoton+MET: prompt
If there is more than one flavor of χ i (N χ > 1), a heavier χ h can decay to the lightest χ l and a photon: χ h → χ l γ. If χ h is pair produced at the collider, diphoton+MET can be important in probing the multi-flavor scenario. For this channel, we follow the diphoton search at 8 TeV of energy and 20.3 fb −1 of data performed by the ATLAS collaboration [75]. The search is aimed to constrain the parameter space of gauge mediated supersymmetry breaking models. The event selection applied in this analysis is as follows. At least two photons with p γ T > 75 GeV and |η γ | < 2.37 are required. Jets are reconstructed using the anti-k t algorithm with a radius parameter of 0.4. The jets are required to have p j T > 30 GeV and |η j | < 2.8. No vetoes are applied on the number of leptons and jets. Several new variables are introduced as follows. An angular separation variable, φ min γ , is defined as the minimum azimuthal angle between E miss T and the two selected photons. In presence of jets, a variable φ min jet is introduced and defined as the minimum azimuthal angle between E miss T and the two highest reconstructed jets, where the jets are required to have p j T > 75 GeV. The total visible energy, H T , is defined as the scalar sum of the transverse momenta of photons, jets and leptons.
There are three relevant signal regions defined in the first three columns of Table 3.
The region WP1 is preferred for a larger mass splitting between m χ h and m χ l , because the cut H T > 400 GeV can only be satisfied with high p T photons. In this signal region the SM background is smaller than in the other two regions, therefore a stricter bound can be achieved. For a lower mass splitting, the region MIS is preferred because it does not have any requirement on the H T value. At 8 TeV, the region WP2 provides a weaker bound compared with the other regions, because of an upward fluctuation in the number of observed events. Note that none of the signal regions veto on the leptons and the jets. In the benchmark model used in the ATLAS analysis, the NLSP are produced in the decay chain of either gluino or chargino, and either jets or leptons are always present in the final states. At the partonic level, the final states of our model contain only a pair of photons and missing energy from the DM particles χ l . Hence vetoing jets and leptons increases the sensitivity of the search to our model.
The main backgrounds for these regions are W γγ, Zγγ and "QCD". In this case, "QCD" is defined as the sum of multi-jet, γ + jets and γγ + jets processes. When estimating the background, we multiplied the LO cross section of the W γγ (Zγγ) background obtained from MadGraph5 with a K-factor of 8.1 (1.8) obtained from VBFNLO [68]. This result still needs to be multiplied by a factor of 3.7 (1.1) to match the distribution reported by ATLAS. 8 The sum of the W γγ and Zγγ missing energy distributions in the WP2 signal region is shown in Fig. 8. In this region the W γγ background dominates Zγγ by an order of magnitude.
For the 14 TeV simulation, we use cuts inspired by the MIS region of the 8 TeV ATLAS analysis, see Table 3. The signal region is defined to have a slightly tighter missing tranverse energy than the 8 TeV analysis. Additionally, in order to increase the sensitivity of the search, jet and lepton vetoes are applied. In order to take into account the low p T jet veto correctly, we generated jet matched samples for both the signal and backgrounds. Contrarily to the 8 TeV case, at 14 TeV Zγγ is the leading background, as a consequence of the lepton veto: the expected background for the 14 TeV signal region is of 2.64 events from Zγγ and 0.26 events from W γγ, assuming 300 fb −1 of luminosity. In the 8 TeV ATLAS analysis, the Zγγ background was estimated using MC, with an associated systematic uncertainty of 50%. Similarly to the previous subsections, we assume that the systematics will improve with luminosity as L −1/2 . The systematic error is thus estimated to be of 0.69 events at 300 fb −1 and 2.2 events at 3 ab −1 . We neglect the systematics for the subdominant W γγ background, since the number of events expected from W γγ is smaller than the uncertainty on Zγγ. Finally we did not simulate the "QCD" backgrounds, but they are expected to be small due to the jet and lepton vetoes together with a harder MET cut.  Figure 8. Comparison between the ATLAS result and our MC simulation for the sum of W γγ and Zγγ in the signal region WP2 (described in Table 3).

LHC diphoton+MET: displaced
For the search of the displaced photon signal plus missing energy, we follow the nonpointing photon analysis in [77], performed by the ATLAS collaboration on about 20 fb −1 of 8 TeV data. The full search also uses the delayed photon measurement, however, due to the complication of modeling the time of flight of the photon from the displaced vertex to the electromagnetic calorimeter (ECAL), we only focus on the measurement of ∆z γ of nonpointing photons (see Fig. 9). For DM signals given by the long-lived χ h → χ l γ decay, ∆z γ can be related to the χ h decay length d in the lab frame: wherer T,z represent the transverse and longitudinal components of the unit vectorr, respectively, as shown in Fig. 9. To obtain the ∆z γ distribution of the DM decay, we first simulate the prompt process, p p → χ hχh , χ h → χ l γ,χ h →χ l γ in MadGraph5, apply the cuts performed in the ATLAS analysis, and reweight the events using the dark penguin form factors. Then we calculate the proper lifetime of χ h and boost it to the lab frame using the momenta of each parton-level event. The angular information of the photon and χ h allow us to calculate ∆z γ in Eq. (3.4) as a function of the decay length. Using this, each simulated MC event contributes to the differential cross section in ∆z γ as where the µ characterizing the probability distribution dP/d∆z γ of the decay is defined as Summing the distributions derived from all the simulated events we obtain the differential cross section in ∆z γ , see Fig. 9. The ATLAS search requires at least two loose photons with |η| < 2.37 and E T > 50 GeV. At least one photon is required to be in the barrel region |η| < 1.37. To avoid collisions due to satellite bunches, both photons are required to have an arrival time at the ECAL t γ smaller than 4 ns, with zero defined as the expected time of arrival for a prompt photon from the primary vertex. We do not attempt to fully simulate t γ , which would require a more complex detector description, but rather we approximate t γ with the time of flight of the χ h , requiring it to be smaller than 4 ns. In our estimation we do not include the detailed isolation cuts on the photon. We also neglect the effect of the displaced decay on the angular acceptance of the photons, simply imposing the requirements on |η| at the level of the prompt event. The signal region also requires / E T > 75 GeV. Finally, to simplify the discussion we assume that every event has a reconstructed primary vertex in the geometrical center of the detector.
For events where only one photon satisfies |η| < 1.37 (i.e. it is in the barrel calorimeter), this photon is used for the measurement of ∆z γ . For events where both photons are in the barrel, the photon with larger t γ is used. We approximate this timing condition by taking the photon emitted by the more boosted χ h , in which case the average decay is more delayed. In Fig. 9 the generated ∆z γ signal distribution is shown, on top of the expected background. The latter is taken from Fig. 4 of the ATLAS paper [77]. Because we are focusing on the non-pointing photon signals, to set constraints on the DM couplings we remove events with |∆z γ | < 30 mm. In our exploratory analysis we only consider the statistical uncertainty on the background, neglecting the effect of systematics.

100 TeV collider and ILC
To give an idea of how much higher energy colliders can improve the sensitivity to the dark penguin, we estimate the bound from the mono-jet search at a 100 TeV collider. Following the discussion in [78], we simulate the DM signal and backgrounds with √ s = 100 TeV and impose the cuts p T > 2.5 TeV, |η| < 2.2 on the leading jet, together with / E T > 3 TeV. A second jet with p T > 100 GeV is allowed as long as it has |η| < 4.5 and the azimuthal separation from the leading jet is ∆φ < 2.2. Events with leptons (taus) are vetoed if the lepton satisfies |η| < 2.5 and p T > 20 (40) GeV. We assume a 2% systematic uncertainty when calculating bounds with 3 ab −1 of data.
To show the importance of lowering the SM background for the dark penguin search, we also estimate the monophoton bound from the ILC-500P and ILC-1000P scenarios in Sec. 4. Our analysis follows the one in [79] by assuming a 500 GeV (1 TeV) ILC with 250 (500) fb −1 of data, and polarizations equal to P − = + 0.8 and P + = 0.5 for the electron and positron beams, respectively. The Z-related SM background e + e − → Z(νν)γ can be eliminated by cutting away the Z pole, but the background process involving a t-channel W cannot be reduced by simple kinematic cuts. The polarization of the ILC beams plays an important role in reducing the latter background, since the W only couples to the left-handed electron. The proposed search requires E γ > 8 GeV, |cos θ γ | < 0.995, and imposes a veto on events with photon energy 238 < E γ < 245 (490 < E γ < 495) GeV at ILC-500P (-1000P) in order to suppress the Z(νν)γ background. Notice that in our simulation of the dark penguin monophoton signals at the ILC, we only considered the case where the photon is emitted by ISR. However, in principle the contribution from the Rayleigh operatorχχB µν B µν is of the same order and should be consistently included. This could lead to quantitative changes in our results, although the qualitative features would remain the same. We leave the inclusion of the Rayleigh operators for future work.

Collider phenomenology of the N χ = 1 case
In this section we present the collider reach on the scenario where only a single species of dark fermion χ is accessible at colliders. In this case, signals in the monojet and monophoton plus missing energy channels can be produced through initial state radiation. In addition, when the DM mass is less than half of the Z mass, the measurement of the invisible Z decay at the Z-pole also sets a strong constraint.

Bound from the Z invisible width
If the DM is lighter than m Z /2, the decay Z →χχ is kinematically allowed and constrained by the measurement at LEP/SLD of the invisible decay width of the Z. For m Z M , the decay Z →χχ is well described by an effective dipole interactionχσ µν χB µν (see Eq. (A.12)), leading to the following approximate expression for the width The uncertainty on Γ(Z → invisible) is of 1.5 MeV [80], therefore the 95% bound reads This bound is shown as an orange curve in the right panel of Fig. 11.

LHC and future colliders
At the LHC, the most promising search for the dark penguin with N χ = 1 is monojet, in which a high p T jet is produced through the QCD ISR process, see TeV run can only cover part of the perturbative parameter space for relatively large Y 4, in which case the perturbative bound is λ √ Y N ∼ < 8π. To compare the dark penguin to the EFT, we also compute the bounds by parameterizing the DM-SM coupling with the effective dipole interaction. The resulting constraint, shown by the dashed blue lines in Fig. 11, is weaker than the dark penguin one, since the dipole EFT neglects the enhancement of the form factor F σ for q 2 2M , as well as the sizable contribution from the q 2 γ µ term in Eq. (2.2) when q 2 is large (in the EFT the q 2 γ µ piece would correspond to a dimension-6 operator, see the first term in Eq. (A.12)). Furthermore, the optimal MET cut for the dark penguin signal is typically softer than for the EFT coupling. In Fig. 12 we show a comparison of the best MET cut for the dark penguin and the effective dipole coupling, as a function of the mediator mass. When the mediators are light, the best MET cut for the dark penguin is lower than for the effective coupling. Notice that in both cases the optimal MET cut increases at larger luminosity, because in our data-driven analysis the systematic uncertainties are mainly determined by the statistics of the control samples.
The LHC monophoton search gives a much weaker bound compared to monojet, and will not be discussed in detail.
Given the limited sensitivity of the LHC to the N χ = 1 scenario, it is important to estimate how the bounds could be improved at future colliders. The study of DM production at very high energy machines, such as a 100 TeV pp collider, mandates the full inclusion of the loop form factors, since the typical partonic energy is much larger than the masses of the mediators. Differently than in the search for heavy particles, however, the significance of the dark penguin signal will not be appreciably improved at a higher energy machine, unless the systematic uncertainty can be greatly reduced. In fact, when the mediators have a sub-TeV χ Z, γχ χ Z, γχ  In v is ib le Z d e c a y Figure 11. Upper bounds on the DM-mediators coupling for N χ = 1, as a function of the mediator mass M . The dark matter mass is assumed to be m χ = 50 GeV (left) and 10 GeV (right). Here we optimize the LHC bound by choosing the best MET cut for each mediator mass. The blue dashed curve shows the bound from the effective dipole interaction,χσ µν χB µν in Eq. (A.12). The monophoton constraint from the ILC500-P (ILC1000-P) assumes 250 (500) fb −1 of data with 500 GeV (1 TeV) center of mass energy and a polarization P − = + 0.8 and P + = 0.5. The orange curve shows the constraint from the current invisible Z decay measurement. The perturbative bound on λ √ Y N depends on the hypercharge of the mediator, λ mass, both the SM background and DM production at 100 TeV are mediated by effectively nearly-massless particles. This results in similar shapes of the MET distribution for the signal and background, and as long as the background is systematics-dominated, the ratio between the signal and background does not vary significantly when increasing the collider energy. The pink curves in Fig. 11 give an estimate of the 100 TeV reach (the search is described in detail in Sec. 3.6). The result is an improvement of the bound on λ √ Y N by only a factor ∼ 2 with respect to the projected 14 TeV LHC constraint. The sensitivity of the 100 TeV collider ameliorates if the MET cut can be kept very weak, below the TeV scale, as in this case the signal region includes the enhancement of the form factors for q 2 ∼ 2M , thus increasing the signal yield. However, it is not clear whether such a low MET cut would be achievable, and if so, whether the systematic uncertainties would be increased beyond the 2% that we are  Figure 12. The MET cut of the projected 14 TeV search that optimizes the DM signal significance as a function of luminosity. As can be seen in the plots, dark penguins with light mediators prefer lower MET cuts. When the mediators are heavier the MET distribution becomes harder, similar to the EFT case, and a tighter cut is preferred. The optimal MET cut increases with luminosity, because in our data-driven analysis the systematic uncertainties are mainly controlled by the statistics of the control samples.
using as benchmark in our estimate. More detailed studies are necessary to precisely assess the reach of a 100 TeV machine.
In comparison to hadron colliders, the monophoton search at the ILC appears to be more promising for the dark penguin signals. As discussed in Sec. 3.6, thanks to the beam polarization and full reconstruction of the missing invariant mass, the SM background can be efficiently reduced. As is shown in Fig. 11, the ILC500-P with 500 GeV center of mass energy, 250 fb −1 of data and polarized beams can improve greatly the sensitivity to the dark penguin.
5 Collider phenomenology of the N χ = 2 case with anarchic DM couplings Differently from the direct and indirect detection experiments, colliders are capable of probing extended dark sectors, through the production of the additional states that accompany the DM. In this section we study the monophoton, diphoton, and rare Z decay constraints on the N χ = 2 dark penguin model, where one heavy dark fermion χ h is present in addition to the light DM particle χ l . In this scenario, large photon signals can be generated through the pair production of dark fermions, followed by the decay χ h → χ l γ. Under the assumption of anarchic couplings between the dark fermions and the mediators, we can compare the bounds from different searches and provide complementary information to the direct and indirect detection experiments. In the diphoton search, the kinematic variable M T 2 can be used to determine the mass of χ h .
The constraint from the invisible Z decay width applies also in the N χ = 2 case, as long as m χ l < m Z /2, as we assume. In addition, when also m χ h < m Z /2, the search at LEP1 for events containing two photons plus missing energy sets a much stronger constraint on the DM couplings. Therefore we focus first on the case m χ h > m Z /2, and discuss separately the more constrained possibility of a lighter χ h , in Sec. 5.3.

LHC monophoton
The monophoton signal arises through the production of χ h χ l , followed by the decay χ h → χ l γ, see Fig. 13. The branching ratio for the decay is either ∼ c 2 w 0.8 or unity, depending on whether the χ l Z channel is kinematically open or closed. Therefore the cross section is much larger than for the ISR process discussed in the case N χ = 1, whose amplitude is comparatively suppressed by the three-body phase space and by the additional electromagnetic coupling. In this section we calculate the bound on λ √ Y N from the 8 and 14 TeV monophoton searches, and study the dependence of the bound on the mass splitting between χ h and χ l . The result is shown in the left panel of Fig. 14. For light mediators, the slope of the bound on λ √ Y N decreases. This can be understood by noticing that the requirement of a minimum photon transverse momentum p γ T imposes q 2 2p γ T , while the dark penguin peaks at q 2 ∼ 2M . Therefore, for M p γ T a significant fraction of the signal does not pass the event selection. The effect is more evident in the 14 TeV case, where the cut on the photon p T is stronger (300 GeV).
To gain some insight about how the mass difference between χ h and χ l affects the monophoton signal, we fix the DM mass to m χ l = 10 GeV and show the photon p T distribution for various choices of m χ h in the right plot of Fig. 14. The coupling λ √ Y N is fixed to the value corresponding approximately to the reach of the 14 TeV LHC. We see that the signal yield decreases significantly only for m χ h 20 GeV. Thus the monophoton search has sensitivity even for fairly degenerate dark fermions, with mass splitting m χ h − m χ l much smaller than the required photon p T , thanks to the boost of χ h from the production process.
In Fig. 15 we show a comparison of the bound obtained using the dark penguin and the one computed using the effective dipole couplingχ h σ µν χ l B µν + h.c. (see Eq. (A.12)). For each mediator mass we vary the photon p T cut between 300 and 1000 GeV in 100 GeV steps, and choose the value that gives the strongest constraint (therefore the dark penguin bound shown in Fig. 15 is slightly better than the one in Fig. 14, which was obtained with a fixed cut p γ T > 300 GeV for all mediator masses). Similarly to the monojet case, the bounds computed by taking into account the full dark penguin are stronger than the EFT results, due to the form factor enhancement at q 2 ∼ 2M and the contribution of the q 2 γ µ term in Eq. (2.2) for large q 2 .  Figure 15. A comparison between the monophoton bounds computed using the dark penguin amplitude (blue) and the effective dipole interactionχ h σ µν χ l B µν + h.c. (orange). For each mediator mass, the bounds are optimized by choosing the best photon p T cut between 300 and 1000 GeV.

LHC diphoton+MET
If the χ h is light enough to be pair produced with sizable rate, the diphoton+MET signal can be generated through the process shown in Fig. 16. Differently from the monophoton case, there can be a partial cancellation of the transverse momenta of the two χ l particles, leading to a suppression of the MET. Nevertheless, the low SM background makes diphoton+MET a promising search to look for signals of the dark penguin. 9 9 If kinematically allowed, the outgoing photons in Fig. 16 can also be replaced by Z bosons, although the smaller BR(χ h → χ l Z) ∼ s 2 w and the additional branching ratios for the Z decays further suppress the cross The three relevant signal regions for 8 TeV are defined in the first three columns of Table 3. For each point in parameter space we compute the constraint from each of the three regions, and choose the strongest one. The reach of the 14 TeV LHC is estimated using the signal region defined in the last column of Table 3. The results for both 8 and 14 TeV are shown in Fig. 17a.
If a signal is observed, the diphoton+MET channel offers the possibility to measure m χ h through the stransverse mass variable M T 2 [81][82][83]. This variable is constructed from the two photons and the missing energy and is defined as where In the above equations µ N is an unknown trial mass, / p T is the transverse missing momentum, p γ j T is the transverse momentum of photon γ j , while the two transverse energies are defined as If the trial mass µ N is chosen equal to the χ l mass, the distribution of M T 2 has an edge at the value of m χ h . In Fig. 17b we show the M T 2 distribution at the 14 TeV LHC with 3 ab −1 of data, for some illustrative choices of m χ h , assuming µ N = m χ l . For each m χ h , we chose the value of the coupling λ √ Y N such that a 3σ excess would be observed at an earlier stage of LHC running, namely with 100 fb −1 of data. The edge at M T 2 ∼ m χ h can be seen in all cases, with some uncertainty of O(10) GeV due to detector effects.
In Fig. 17b we assumed the mass of the DM to be known, which allowed us to set µ N = m χ l in the computation of M T 2 . But how could m χ l be determined experimentally? References [84,85] showed that the value of m χ l can be estimated by observing a kink of the edge of the M T 2 distribution as a function of µ N . However, a large number of signal events is required in order for the kink to be observable, making the application of this method, and thus the determination of the DM mass directly in the diphoton process, not likely for our model. On the other hand, direct detection constraints may hint that the DM is light. As we will discuss in Sec. 7, in the region of parameters where the LHC has sensitivity to a dark penguin signal, current direct detection searches require m χ l 10 GeV. Hence one can argue that if a diphoton+MET signal is observed at the LHC, and interpreted as involving as missing energy the particle that provides the dominant dark matter density in the universe, one should use the value µ N 10 GeV.
It is interesting to compare the results from the monophoton and diphoton searches. The comparison requires some further assumption, since the monophoton search sets a constraint on the dark penguin coupling between χ h and χ l , whereas the diphoton search constrains the coupling between χ h and χ h . As discussed in Sec. 2, here we focus on the totally anarchic scenario, where all couplings are assumed to be of the same order. In Fig. 18 we show a comparison between the 14 TeV monophoton and diphoton searches, for the illustrative choice (m χ l , m χ h ) = (10, 300) GeV. Diphoton gives a slightly stronger bound, but the difference is below the uncertainties associated with our analysis. Notice that, while it is useful to compare the two channels within a specific natural scenario such as the anarchic one, in general the monophoton and diphoton searches are complementary to each other.

Constraint from rare Z decays
We now turn to the case of light χ h . If m χ h < m Z /2, the decay Z → χ hχh → γγχ lχl is kinematically allowed, and is constrained by LEP1 data collected at the Z-pole. A search performed by OPAL [86] for events containing at least two photons and missing energy sets a strong constraint on this decay. The search region required exactly two photons, satisfying m Χ l 10 GeV the cuts E γ > 1 GeV, |η γ | < 1.74, m γγ > 10 GeV and π − ∆φ γγ > 0.0873 (acoplanarity angle). In addition, at least one of the two photons was required to have |η γ | < 1.10, and additional photons with E γ > 0.5 GeV were vetoed. No events of this type were observed in 43 pb −1 of data, corresponding to 1.8 × 10 6 Z bosons. The expected SM background comes from e + e − → Z(νν)γγ, where the two photons arise from initial state radiation, and is quoted by OPAL to amount to 0.2 events. 10 Therefore, assuming no systematic uncertainty, the 95% CL limit is of 2.8 signal events. In order to set a bound on the parameter space, we still need to compute the efficiency of the cuts on the signal. For m χ h = 30 GeV the efficiency is of ∼ 70%, leading to the 95% CL bound (the branching ratio for the decay χ h → χ l γ is unity, which under the assumption of anarchic couplings is much stronger than the other LEP constraint from the invisible Z width. The corresponding excluded region is shown in Fig. 19, where it is compared to the projected sensitivity of monophoton and diphoton searches at the 14 TeV LHC. We conclude that if m χ h < m Z /2, only the diphoton search will be able to marginally improve the LEP constraint. This result motivates our focus on the scenario m χ h > m Z /2, in which case LEP data only give a mild constraint from the invisible Z width. m Χ l 10 GeV 6 Non-pointing photon signals As discussed in Sec. 2.1, when the flavor structure of the dipole operator is aligned with the mass matrix, the decay rate of a heavier fermion χ h into the DM χ l and a photon can be highly suppressed, while the coupling of γ, Z to a pair of χ h 's or χ l 's is still sizable. This motivates the study of displaced photon signals from pair produced χ h 's, each decaying into χ l + γ. The bounds obtained using the ATLAS non-pointing photon search [77] are shown in Fig. 20, for some representative choices of the parameters. The bounds are computed by means of a simple counting experiment, by comparing the signal and background yields in two different regions, corresponding to |∆z γ | > 30 mm (exclusion shaded in blue) and |∆z γ | > 220 mm (shaded in red). The displacement ∆z γ was defined in Fig. 9. For both signal regions, an upper bound |∆z γ | < 750 mm is also imposed, since the ATLAS paper does not report the background expectation for larger values of the displacement. In the first signal region, approximately 149 events were expected from the SM background, and 140 were observed. In the second region, approximately 14 events were expected and 18 were observed. A few comments are in order. The first observation is that, independently of the value of m χ h , the lower bound on the coupling 11 λ as a function of scales as λ ∝ −1/4 . This can be understood by noticing from Eq. (3.5) that, for large Γ χ h (the lower bound on λ corresponds to large displacement of the decay), the number of signal events passing the cuts is N S ∝ λ 4 Γ χ h , since the production cross section scales like λ 4 . Recalling that Γ χ h ∝ λ 4 2 and that in our counting experiment the exclusion bound is given by N S = constant, we obtain the scaling λ ∝ −1/4 . Second, comparing the cases m χ h = 300 GeV and m χ h = 50 GeV we note that the bound is weaker in the latter case, especially for the harder cut |∆z γ | > 220 mm. The reason for this is twofold. In first place, the efficiency of the cuts on the E T of the photons and on the / E T (see Sec. 3.5 for details) is smaller for lighter χ h . In second place, a lighter χ h is more boosted, which implies that the photon from its decay is typically more collinear with χ h . This fact, combined with the upper cut of 4 ns that we impose on the time of flight of χ h to avoid spurious collisions due to LHC satellite bunches, implies that the typical photon displacement is smaller for lighter χ h . As a consequence, for m χ h = 50 GeV the cut |∆z γ | > 220 mm significantly reduces the signal yield.

Comparison to the direct and indirect detection results
The comparison between the collider and direct detection constraints is shown in Fig. 21, under the assumption of anarchic couplings between the mediators and the dark fermions.
Here we follow the analysis in [61,87], comparing the bounds set by the two classes of experiments on the magnetic dipole moment µ χ , defined by the effective operator where χ is the DM, and related to the parameters of the dark penguin that couples two DM particles to the photon by Since the momentum exchange in direct detection experiments is small, we can safely integrate out the mediators and use the effective description in Eq. (7.1). Results from direct detection experiments were obtained in [61], with the assumption of the standard halo model and velocities v esc = 544 km/s, v = 232 km/s, v 0 = 220 km/s. For the local DM density, the standard value ρ = 0.3 GeV/c 2 /cm 3 was assumed. Figure 21 shows the 90% CL bounds from CDMSlite [88], SuperCDMS [89], the XENON10 S 2 -only analysis [90], XENON100 [91], LUX [92], CDMS-II-Si [93], CDMS-II-Ge low threshold [93], and CoGeNT2014 data [94], together with the 68% and 90% CL allowed regions for DAMA [95] (assuming quenching factor Q Na = 0.30), CoGeNT2014, and CDMS-II-Si. For XENON10 we take the result with a conservative setting of the electron yield to zero below 1.4 keVnr, as in Ref. [91]. For LUX we adapt the limit corresponding to zero observed events. For further details on the experimental data, as well as on the assumptions made on the low energy thresholds and quenching factors, we refer the reader to Ref. [61]. We discuss first the N χ = 2 case, shown in the top row of Fig. 21. The red horizontal line indicates the upper bound on µ χ obtained from the invisible Z decay at LEP, whereas the green and cyan horizontal lines correspond to the bounds obtained from the projected monophoton and diphoton searches at the 14 TeV LHC, respectively. Due to the anarchic assumption, all these searches can be interpreted as effectively constraining the size of the dark penguin that couples two DM particles to the photon, and thus the size of the effective interaction in Eq. (7.1) via Eq. (7.2). This allows us to make the comparison to direct detection experiments.
To gain some understanding of the dependence of the bounds on the mass of χ h , we show results using two benchmark values, 50 and 300 GeV. On the other hand, the collider bounds are essentially independent of m χ l in the range considered m χ l 20 GeV, the only appreciable effect being that the LHC monophoton and diphoton bounds worsen slightly for m χ l 10 GeV in the case m χ h = 50 GeV. The plots in Fig. 21 assume the benchmark value of the mediator mass M = 500 GeV, but they serve as useful order of magnitude estimates also for heavier mediators. By the end of the high-luminosity LHC run, monophoton and diphoton+MET searches will be able to test values of the dipole interaction strength that are comparable to the current XENON10 reach for m χ l > 10 GeV, and will provide a much better bound than the current direct detection limits for lighter χ l .
The N χ = 1 case is shown in the bottom panel of Fig. 21. The projected 14 TeV LHC monojet bound (purple) is slightly worse than the current limit from the invisible Z width (red), showing that to improve the constraint it will be necessary to go to future colliders. The 100 TeV collider monojet bound is shown in cyan and corresponds to µ χ ∼ < 2 × 10 −17 e cm, while the ILC1000-P monophoton search (green) can reach a bound 5 × 10 −18 e cm. For comparison, the LEP constraint derived in [42] corresponds to µ χ ∼ < 10 −16 e cm. Similarly to the N χ = 2 case, collider searches play a more important role than direct detection experiments when the DM is light, m χ ∼ < 10 GeV. The collider bounds can also be compared to the indirect detection searches, where a pair of DM particles annihilates into photons. There are two mechanism through which the    Figure 21. Top row: comparison between the projected collider constraints and the current bounds from direct detection experiments for N χ = 2, based on the assumption of anarchic couplings between the mediators and dark fermions. Here we adapt the direct detection bounds from [61,87], assuming a standard halo model (SHM). See the text for more details. The magnetic dipole moment is defined by the effective coupling of the DM to the photon, (µ χ /2)χσ µν χF µν . The red solid lines show the bound from the measurement of the invisible Z decay width at LEP. The green (cyan) lines show bounds from the monophoton (diphoton) searches at the 14 TeV LHC, with solid (dashed) line for 300 (3000) fb −1 of data. Bottom: projected collider constraints and current bounds from direct detection experiments for N χ = 1. The lines correspond to the invisible Z width at LEP (red), 14 TeV LHC monojet (purple), 100 TeV collider monojet (cyan) and ILC1000-P monophoton (green). annihilation can proceed in our model: either via two dipole interactions connected by a t-channel DM exchange, or via a Rayleigh operator. Using naive dimensional analysis, the ratio of the former amplitude to the latter [41] is ∼ λ 2 M/(16π 2 m χ ). Considering the most conservative choice of parameters relevant to our discussion, namely λ ∼ 4, M ∼ 200 GeV and m χ ∼ 10 GeV, we find the ratio to be about 2. We conclude that for our purposes it is a safe approximation to neglect the contribution of the Rayleigh operator, which is typically strongly subdominant compared to the diagram with t-channel DM exchange. The latter  Figure 22. Upper bounds on the annihilation of DM into photons, derived from the projected collider constraints computed in this paper. For N χ = 2 (left) we show the projected LHC bounds, while for N χ = 1 (right) the reach of future colliders is presented. In the right panel, the Fermi bound corresponds to the expected limit from the 3.7 years data, as reported in [96] (we take the most stringent limit in Ref. [96], obtained using the R16 Einasto DM profile).
dominates the annihilation, with a cross section [40] σχ χ→γγ v = 1 2π µ 4 χ m 2 χ . (7.3) As shown in Fig. 22, if more than one flavor of dark fermions are accessible at the LHC, the size of the annihilation into photons can be tested down to approximately 10 −33 cm 3 /s. In the N χ = 1 case, the ILC1000-P can set the strongest bound on the annihilation, of order 10 −32 cm 3 /s.

Conclusions
The search for DM plays a central role in the physics program of current and future collider experiments. For many plausible DM-SM couplings, a correct description of the DM production processes, which is crucial to extract meaningful information from the experimental results, requires that the EFT parameterization be UV completed by a simplified model with light mediators. This is especially important when the DM-SM coupling arises at loop level. In this paper we performed the first simplified model collider study for a loop process, focusing on the Dark Penguin, whose form factors reduce to the magnetic dipole operator at low energies. We computed bounds from monojet, monophoton, and diphoton searches at the 8 and 14 TeV LHC, as well as from the future ILC and 100 TeV hadron collider. Differently from searches for EFT interactions, when light mediators are included the optimal search strategy requires the cuts on the MET to be as loose as possible, to capture the enhancement of the form factors near the threshold for production of on-shell mediators. As we showed through a detailed comparison, for light DM mass the collider bounds are complementary to those derived from direct and indirect detection experiments.
Based on general considerations, it is plausible that the dark sector may be endowed with a non-trivial flavor structure. If some of the additional states beyond the DM are kinematically accessible, collider experiments offer a unique opportunity to probe the dark flavor. By employing a simplified model with a second dark fermion in addition to the DM, we showed that, under the natural assumption of anarchic structure in the dark sector, the bounds on the DM-SM coupling set by collider searches are much stronger than in the case where the DM is the only accessible state. Collider searches not only have the capability to probe the flavor structure of the dark sector, but might even allow to measure some of its properties. For example, the determination of the mass of the heavy dark fermion could be possible by using M T 2 in the diphoton+MET channel. Furthermore, if the flavor structures of the dark penguin and the dark fermion mass matrix are nearly aligned, the decay of the heavy dark fermion into DM and photon can be displaced. In this case, it is achievable to extract information about the small mixing angle in the dark sector from the search for non-pointing photons+MET.
We end with an outlook to future developments. While our study was focused on the dark penguin mediating the dipole operator, our method is fully general and can be extended to any other loop-mediated interaction between the DM and the SM fields, such as for example the Rayleigh andχχG µν G µν operators, where G µν is the SM gluon field strength. It would also be interesting to analyze in detail the prospects for the ILC, which the preliminary results presented here show to be very promising. Additionally, while our first estimates for a very high energy hadron collider appear less favorable, the detailed assessment of the design requirements that would allow a substantial improvement warrants a dedicated study. and we defined d 3 x ≡ 1 0 dxdydz δ(x + y + z − 1). It is straightforward to check numerically that the Eqs. (A.9) agree exactly with Eqs. (A.1) for q 2 < 4M 2 , whereas for q 2 > 4M 2 the Feynman parameter integrals can become numerically unstable, and only the expression in terms of scalar integrals should be used. 13 By using Eqs. (A.9), one can derive the expressions of the form factors for large mediator mass, M f = M s = M m i , m j , q 2 , In this limit, the amplitude in Eq. (2.2) can be seen as generated by the effective Lagrangian

B Statistics
In this appendix we briefly describe our procedure for setting limits on the DM parameter space. Given the expected number of background events N B , we exclude a signal model yielding N S events if where dP/dx (x; N B + N S ) is the normalized probability distribution function for signal plus background, N obs is the number of events observed, and p is the chosen probability. For example, for a 95% CL exclusion, p = 0.05. Notice that we are setting a 'one-sided' limit. Throughout our analysis we neglect systematic uncertainties on the signal, because they are subleading to those associated with the background. Under the assumption that signal plus background follows a Poisson distribution with mean N B + N S , and neglecting all systematic uncertainties, the exclusion limit is given by Γ(N obs + 1, N B + N S ) Γ(N obs + 1, 0) = p , (B .2) where Γ(s, q) = ∞ q t s−1 e −t dt . For example, for p = 0.05, we find for N B = N obs = 0 that N S 3.0.
If we assume that signal plus background follows a Gaussian distribution with mean N B + N S and standard deviation N S + (δN B ) 2 , where δN B = N B + ( N B ) 2 , with the relative systematic uncertainty on the background, we find the exclusion limit 3) 13 We have cross-checked our results against Ref. [49], where the form factors were computed for the special case mi = mj, Y = 1/2, N = 2. We find agreement, except for a few small differences that we believe are due to typos in App. A of [49]: 1) in their Eq. (A7), a factor z should multiply the second term on the right-hand side (RHS), and an overall factor (−1) should multiply the RHS; 2) in their Eq. (A8), an overall factor (−2) should multiply the RHS.