UvA-DARE (Digital Academic Repository) Effective field theory of dark matter: a global analysis

: We present global (cid:12)ts of an e(cid:11)ective (cid:12)eld theory description of real, and complex scalar dark matter candidates. We simultaneously take into account all possible dimension 6 operators consisting of dark matter bilinears and gauge invariant combinations of quark and gluon (cid:12)elds. We derive constraints on the free model parameters for both the real ((cid:12)ve parameters) and complex (seven) scalar dark matter models obtained by combining Planck data on the cosmic microwave background, in principle explain the Galactic centre excess, but that such interpretation is in strong tension with the non-detection of gamma rays from dwarf galaxies in the real scalar case. In the complex scalar case there is enough freedom to relieve the tension.


Introduction
Overwhelming observational evidence points to the existence of dark matter (DM) , from Galactic up to cosmological scales [1][2][3][4] The most widely discussed and well-motivated particle DM candidates are weakly interacting massive particles (WIMP), which arise from some of the most popular extensions of the SM, and naturally lead to the right DM abundance [4].
Well-motivated ultraviolet-complete particle-physics models accommodate good DM particle candidates in their mass spectrum. This is the case for example of supersymmetric theory with R-parity conservation [2] or Universal Extra Dimensions theories [5]. On the other hand, one can adopt a model independent approach that makes minimal assumptions on the DM particle and its couplings with SM particles. In the framework of Effective Field Theories (EFT) for particle DM, DM would be the only additional degree of freedom beyond the SM accessible by current experiments [6]. Therefore, the interactions of the DM particle with SM particles are described by effective operators (of dimension 6 or higher). Those can be predictive if the energy scale of the experiment under investigation is lower than the energy scale of the operator's coefficients, while will break down once the energy scale of the experiment is of the order of the mass of any particle mediating the DM-SM interactions.
Several detection strategies are applied in order to detect the elusive nature of DM. Direct detection searches look for the recoil energy of nuclei scattered off by DM particles in large, underground laboratories. Indirect detection searches aim to detect the final stable products (like gamma rays, neutrinos or charged cosmic rays) of DM annihilation or decay JHEP09(2016)077 above the large astrophysical background. Finally, searches for DM at colliders are based on the possibility to look for the production of new particles beyond the SM.
EFT have been shown to be able to capture the main features of generic WIMP candidates in their range of validity and being a powerful framework to compare theory against data. EFT operators have been studied extensively in context of various experiments, see for example [7][8][9]. However, in these studies each operator was considered separately to draw phenomenological implications and detectability prospects. In this work, we present the first global analysis of all relevant EFT operators simultaneously in light of the latest constraints from indirect and direct searches for DM. We perform a Bayesian statistical scan of the EFT parameter space, as described in section 2, by implementing all the latest experimental constraints, following an approach similar to that adopted for supersymmetric scenarios in refs. [10][11][12][13][14][15][16][17]. Previously, Balázs et al. [18] presented a similar work with a full set of EFT operators. Our analysis goes beyond that of ref. [18] in several respects: we perform a more thorough statistical analysis that addresses the dependence on priors, and we perform a comparison between profile likelihood and posterior distributions; we include the contribution from DM-gluon operators that was previously neglected while being potentially sizeable; we make use of the latest data from the Fermi Large Area Telescope (Fermi -LAT); and we consider the effect of DM annihilation on the cosmic microwave background (CMB) anisotropies measured by Planck. Blennow et al. [19] also did a global analysis of EFT operators, they studied a Dirac fermion DM candidate while in this work we study a scalar candidate. We do not include collider constraints from the LHC, as these do not rigorously exist for scalar dark matter, and are typically subdominant to direct search constraints for the masses of particular interest to us here [20,21].
This work is organised as follows. In section 2 we describe the theory of EFT operators we use. In section 3 we describe our statistical approach and computational method, in section 4 we present the results, and in section 5 we present our conclusions.

Effective field theory operators
We focus on the interactions of DM with quarks and gluons, since these interactions are currently probed very efficiently by so-called 'direct detection' experiments, and by searches for DM production at the LHC. It is relatively straightforward to extend our formalism to include leptons (e.g. [22]), but we leave such refinements for future work. We further make some well-motivated simplifying assumptions about the nature of the interactions. In particular, we work in the limit in which the particles mediating the interactions between the DM, quarks, and gluons are heavy compared to the energies of interest. In that limit, all theories map into an EFT which encapsulates the interactions of the DM with the Standard Model via a set of non-renormalizable interactions [6-8, 23, 24], where the coefficients λ i have dimensions of inverse mass to the appropriate power such that the over-all dimension of L remains four and the O i are a set of operators consisting of a DM bilinear contracted with a gauge invariant combination of quark and/or gluon fields.

JHEP09(2016)077
Real scalar DM operators Table 1. The EFT operators for real and complex scalar DM interacting with quarks and gluons. Also indicated is the mapping for each coefficient λ i to the notation of [7], which operators contributes to σ SI , and if σ ann v is s-wave or p-wave dominated in the non-relativistic limit.
As a starting point, we focus on the case in which the DM is a single species of (real or complex) scalar particle that is a singlet under the electroweak symmetry, as this limits the number of Lorentz structures describing its interactions. 1 We restrict our basis of operators to those which are leading in the sense of being the least marginal at low energies and which represent the leading structures consistent with the principle of minimal flavor violation (MFV) [25], which dictates that their contributions to flavour-changing neutral currents follows the same CKM structure of the SM itself, mitigating the otherwise extreme constraints from the null searches for non-SM sources of flavour violation. After rotating the quarks into the mass basis, this effectively results in the quark vector bilinears having a generation-independent coupling whereas the scalar bilinears are weighted by the quark mass. The resulting set of operators for real (R) and complex (C) DM are shown in table 1. Also indicated in the table are the operators which make velocity-unsuppressed contributions to spin-independent scattering with nuclei (σ SI ) or annihilation ( σ ann v ) in the non-relativistic limit. Note that the possible annihilation channels for these models are the kinematically available quarks and gluons. The branching ratios are determined by the relative strength of the operators. For each individual operator the ratios to different JHEP09(2016)077 flavours (when applicable) are democratic except for the operators weighted by the quark mass. For these the branching ratios to different quark flavours do have an m 2 q dependence, i.e. the heaviest available is favoured.
While typically an analysis will assume that one operator or another dominates, any realistic UV model of scalar DM will involve several in concert with related coefficients (see [24,[26][27][28][29][30][31][32] for examples). To truly represent the heavy-mediator limit in general, one must allow for combinations of interactions. Combined with the DM mass, this defines a parameter set of five quantities for real DM and seven for complex. This work represents the first truly general analysis of the scalar singlet DM parameter space in the EFT limit.
It is also worth mentioning that the EFT description will fail to accurately describe observables whose typical momentum transfer is large enough to be on the order of the particles mediating the interaction. Particularly for large momentum transfer processes such as at the LHC, this implies that limits derived in an EFT context do not apply to models in which the mediator masses are TeV [33,34].

Statistical framework
We use highly efficient Bayesian methods to explore the models but we present our results both in Bayesian and in frequentist terms. Our approach is based on Bayes' theorem (see e.g. [35]) where D are the data and Θ are the model parameters of interest. Bayes' theorem states that the posterior probability distribution function (pdf) p(Θ|D) for the parameters is obtained from the likelihood function p(D|Θ) ≡ L(Θ) and the prior pdf (or "prior" for short) p(Θ). In this article we are primarily interested in parameter inference, therefore the Bayesian evidence p(D) merely act as a normalisation constant, and will not be considered further in the following analysis. In order to study the constraints on a single parameter of interest θ i , one can consider either the one-dimensional marginal posterior, or the one-dimensional profile likelihood. The marginal posterior is obtained from the full posterior distribution by integrating (marginalising) over the unwanted parameters in the n-dimensional parameter space: On the other hand, the profile likelihood function for θ i , instead, is found by maximising over the parameters that are not of interest: The extension of these concepts to more than one parameter is straightforward. The profile likelihood and the marginal posterior are two different statistical quantities that may lead to different conclusions about the parameter space of interest. The marginal posterior integrates over hidden parameter directions and therefore correctly accounts for volume JHEP09(2016)077 effects; it peaks at the region of highest posterior mass. The profile likelihood peaks at the region of highest likelihood. It is oblivious to volume effects, but is an excellent quantity to find small regions of high likelihood in parameter space. These two quantities do not necessarily lead to the same conclusions for non-Gaussian likelihoods, and the maximum of information about the model parameter space is obtained by studying both of these quantities. Therefore, in the following we present results for both the marginalised Bayesian posterior and the profile likelihood.

Experimental constraints and the likelihood function
The experimental constraints are implemented as a joint likelihood function L with each component representing different contraints.
where L Ωχh 2 is the part corresponding to measurements of the cosmological DM relic density, L DD direct DM detection constraints and L CMB , L dSph , L GCE are from DM indirect detection constraints. We discuss each component in turn: ln L Ωχh 2 : we apply a Gaussian likelihood taken the Planck CMB data constraint on the DM relic abundance. We use as central value the result from Planck temperature and lensing data Ω χ h 2 = 0.1186 ± 0.0031 [36] with a (fixed) theoretical uncertainty, τ = 0.012, to account for the numerical uncertainties entering in the calculation of the relic density.
ln L DD : for DM direct detection we use upper limits from the LUX experiment [37], as implemented in the LUXCalc code [38]. We adopt hadronic matrix elements determined by lattice QCD [39,40]. We use a local DM density of ρ dm = 0.4 GeV cm −3 .
ln L CMB : DM annihilating into ionising particles during the cosmic dark ages will broaden the last scattering surface of the CMB. This would modify the CMB anisotropies as measured by the Planck satellite. As Planck has not measure any such effect [41] we can use the data put limits on our models. We use the likelihood as defined by [42] using the updated analysis from [43].
ln L dSph : as limits from indirect searches for DM, we focus here on limits from gamma ray observations with the Fermi -LAT. As the arguably most robust limits in indirect searches, we adopt constraints that were obtained from the non-observation of a gamma ray signal from a dozen of dwarf spheroidal galaxies. We use here results from [44], which are based on a combined analysis of six years of Fermi -LAT data, and take into account uncertainties in the DM content of each dwarf spheroidal. The results of that analysis were presented as tabulated likelihood functions, which allows us to apply them to models with arbitrary gamma ray spectra.
ln L GCE : given the excitement about a possible gamma ray DM signal from the Galactic centre (see [45][46][47] and references therein), we also include a likelihood that evaluates JHEP09(2016)077 the compatibility between this gamma ray excess and a DM annihilation signal as predicted by our models. To this end, we use the results from [48], which account for systematic correlated uncertainties related to the subtraction of Galactic diffuse foregrounds along the line-of-sight towards the Galactic centre.
Both the dwarf spheroidal and the Galactic centre likelihoods we have used are conveniently packaged in the gamLike 2 code.

Priors
As we are using Bayesian inference we need to provide prior probability distributions for our parameters. For the DM mass parameter m χ we consider the typical range m χ ∈ [1, 1000] GeV advocated for WIMPs, and we adopt a prior uniform on the log of the quantity, which reflecs a state of indifference with regards to the mass scale.
Given a mediator particle with mass M with couplings g 1 , g 2 to the DM and SM particles respectively, we can write the operator coefficients λ i as where k i is an operator specific constant. In this way, conditions on the underlying physics can then translate onto the coefficients. The theory has to be perturbative g 1 g 2 ≤ (4π) 2 , and, for the EFT description to be valid, we, at least, need M > 2m χ . These conditions imply The bound's strong dependence on m χ makes the operator coefficients very hard to scan efficiently. We could define hyperparameters to eliminate the mass dependence, but together with the unavoidable lower bound, this would introduce a prior preference towards higher DM masses. Instead we use the fact that, as we shall see later, the constraint from producing the correct relic density is much more constraining than the theoretical bound. We therefore simply impose λ i < 1. There is no physical motivation for a lower bound, so we conservatively adopt 10 −20 , since if λ i = 10 −20 then the operators O i do not contribute significantly to any observable.
In order to assess the prior dependence of the posterior, we consider here two different kind of priors to scan over λ i . We first consider a prior uniform on the log, i.e. log 10 λ i ∈ [−20, 0]. Then, we define two hyperparameters: a common energy scale to all operators A, and the coefficients f i such that λ i = f i A.

Scanning methodology
To map out the posterior and the likelihood we have developed the EFTBayeS package in which we implemented the Lagrangians of our models in FeynRules [49]. This is then interfaced with a modified version of MicrOMEGAs v2.4 [50] to compute the relic abundance of DM and direct and indirect detection rates. The code LUXCalc v1.0.1 [38] has been used to compute the LUX experiment likelihood and PPPC 4 DM ID [51,52] tables to calculate the photon spectrum of DM annihilation products. We have used tabulated values of f eff from [43]. f eff is the efficiency for DM annihilations to deposit energy into the gas medium in the cosmic dark ages and is used in the CMB constraint.
For the exploration of the EFT models, the EFTBayeS code uses MultiNest v2.18 [53,54] nested sampling algorithm. MultiNest is an extremely efficient scanning algorithm that can reduce the number of likelihood evaluations required for an accurate mapping of the posterior pdf by up to two orders of magnitude with respect to conventional MCMC methods. This Bayesian algorithm, originally designed to compute the model likelihood and to accurately map out the posterior, is also able to reliably evaluate the profile likelihood, given appropriate MultiNest settings, as demonstrated in [12].

Results
In this section we present the impact of aforementioned experimental data on the parameter space of the real and the complex scalar DM candidates. We first show the profile likelihood and posteriors when applying the relic density determination and all current experimental limits. We then show what happens when we also assume that the Galactic centre excess is entirely due to the annihilation of our DM candidates. This will serve as an example to demonstrate the constraining power of an astroparticle detection, when combined with the relic density constraint.

Impact of Planck, LUX, and Fermi -LAT
In figure 1 we show the two 2D marginalised posterior and the 2D profile likelihood for both the real and complex scalar DM candidate in the planes of m χ and either σ ann v or σ SI . This is after including all experimental limits in the likelihood, but excluding the Galactic centre excess. In order to illustrate the dependence of the results on priors, we show the posterior distribution under the Log and the Dirichlet priors as defined in section 3.2.
We begin with discussing the features of in the σ ann v panels. We see that the value of the annihilation cross-section is quite well-constrained, due to the fact that it is directly related to the relic density, which is a well measured quantity. For the complex scalar candidate we have two regions in the posteriors, because the C3 and the C4 operators are p-wave rather than s-wave as the rest of the operators are (all real scalar operators are s-wave.) In the lower region it is one of the p-wave operators that dominates, while in the upper region it is one of the s-wave operators. By a dominating operator we mean that the corresponding coefficient λ i is large enough to provide the correct amount of relic density by itself. The differences we see between the posteriors and the profile likelihood is due to the difference between marginalisation (i.e. integration) and profiling (i.e. maximation). Figure 1. Confidence region and credibility regions for real and complex scalar DM candidates in observable vs mass planes. In the m χ , σ ann v panels we show the 68% contours while in the m χ , σ SI panels we show the 95% contours. This is done for clarity of presentation; the posteriors in certain regions are fairly flat which leads to very noisy contours. The experimental constraints applied are the relic density of DM, limit on spin-independent scattering cross-section from LUX, and annihilation cross-section limits from both the Planck CMB measurement, and the stacked dSph analysis from Fermi -LAT. The different posteriors are using different priors; the Log prior and the Dirichlet prior (α = 0.1) as defined in section 3.2.

JHEP09(2016)077
In the complex scalar case, the values of the annihilation cross-section between the p-wave region and the s-wave region still have a large profile likelihood as there is at least one combination of input parameters for which the operator coefficients give rise to the correct relic density. However, as there are not very many of these possible combinations in the input parameter space, they have a negligible effect on the posteriors.
In the real scalar case the most noticeable difference between the profile likelihood and the posteriors are the downward spikes in σ ann v when the DM mass is close to a quark mass. This is a kinematical effect that enhances the cross-section in the early universe when new annihilation channels open up [55]. In this case the new channels are additional JHEP09(2016)077 quarks in the R1/C1, R2/C2 operators. This effect isn't present in today's colder universe which is why we get the downward spikes in σ ann v today. These are tiny, i.e. highly tuned, regions in the parameter spaces which is why we only see them in the profile likelihood.
The bottom of the two spikes at m χ ≈ 1.3 GeV and 4.7 GeV (i.e. charm and bottom quark masses) avoid all limits. In fact the only term in the likelihood that has any impact for m χ < 5 GeV is the CMB term, whereas the constraining power of the LUX experiment vanishes as the recoil energies of these low mass DM candidates fall below the detection threshold. As for Fermi -LAT dSph results, they should have an impact on these region, but the lowest photon energy (E γ ) considered in the likelihood published by Fermi -LAT and implement here, is 500 MeV, and for m χ 5 GeV the spectra actually peaks below E γ = 500 MeV. For m χ = 5 GeV 75% of the photons has an energy less than 500 MeV, for m χ = 1 GeV that number increases to 99%. 3 We strongly suspect that these low mass regions would be disfavoured in an analysis extended to lower E γ .
This limitation in the dSph likelihood also explains the existence of slightly favoured regions seen in the complex scalar profile likelihood maps for m χ and σ ann v in figure 3 below. There is a slight upward fluctuation in the low energy bins of the dSph likelihood, this does not matter for larger m χ as their spectra is fitted against a large set of bins. However, as stated above, for low m χ much of the spectra is below the lowest energy bin of 500 MeV and the importance of the slight upward fluctuation is enhanced.
In contrast with the annihilation cross-section the spin-independent scattering crosssection is only bounded from above by the relic density. This is because not all operators contribute to σ SI , which means that as long as one of the non-contributing operators' coefficient is large enough to provide the correct relic density, the coefficients of the contributing operators are free to be arbitrarily low. This means that the profile likelihood region actually extends down to arbitrarily low numbers. This freedom to put the operator coefficients at arbitrary low values is also the reason why the two posteriors differ, i.e. the tails of the distributions depend on the choice of priors for the operator coefficients.
This is easier to understand if we look at the 1D posteriors and profile likelihood in figures 2 and 3 for real and complex scalar respectively. In the panels for the operator coefficient λ i we clearly see that the tails are prior dependent. The profile likelihood is in fact flat for low λ i and following Bayes' theorem (3.1), when the likelihood is constant, then the posterior simply traces the prior. We note however that the position of the peaks does not depend on the choice of prior (although the height does, since it depends on the size of the tails through the normalisation of the posterior) and they tell us that the most probable solution is when a single operator provides all the DM relic density by themselves.
For the real scalar candidate low masses are effectively disfavoured by Planck, Fermi -LAT and LUX with the exception of the small regions around 1.3 GeV and 4.7 GeV as discussed earlier. In the complex scalar case the mass is much less constrained due to the presence of p-wave operators.

JHEP09(2016)077
Best fit points for the real scalar DM case  Table 2. Best fit points (i.e. minimal χ 2 ) for both the real and complex scalar DM candidates with and without fitting to the Galactic centre excess. The p-values are calculated only using χ 2 contribution from the Galactic centre excess, under the fairly bold assumption that the test statistic is chi-squared distributed with 24 − 3 = 21 degrees of freedom.

GCE as example of additional measurement
In figures 4 and 5 we show the 1D profile likelihood and posteriors for real and complex scalar DM candidate when we include an actual measurement in the likelihood, in addition to the relic density. The measurement we consider is the Galactic centre excess because it is an actual well-studied measurement, and still compatible with a DM interpretation. Although a more standard astrophysical interpretation, e.g. in terms of millisecond pulsars [56][57][58][59][60][61], is possible or even probable, we assume here that the excess is fully explained by our DM candidates. The impact of including this measurement in the likelihood is drastic. For the real scalar case only the R2 operator is the possible solution with a m χ ≈ 40-60 GeV, whereas the R1 and R3 operators are disfavoured because they contribute to the spin-independent scattering cross-section and are therefore disfavoured by LUX; and the R4 operator favours slightly lower masses which are excluded by the Fermi -LAT dSph analysis.
In fact, even the R2 operator is in tension with the limit from the dSph analysis, but it is the least disfavoured operator. If we look closer at the best fit point from the scan in table 2 we see that it is a good fit to the excess. However, if we compare with the best fit point when not fitting the excess we see that the other data, mainly the dSph, contributes a ∆χ 2 = 5.1 which is sizeable. In summary, if one believes the Galactic centre excess is due to DM annihilation, the tension with the non-detection of gamma rays from dwarf galaxies would disfavour the real scalar model is excluded as the dark matter candidate.
In the case of the complex scalar candidate, the two additional operators, C3 and C4, are p-wave operators and thus cannot explain the excess in themselves. They do, however, relax the connection between the relic density and the annihilation cross-section which somewhat relieves the tension from the dSph limit. If we look at the best fit point for the complex case in table 2 we again see that it is a good fit to the excess. Comparing with the best fit point without the excess we see a ∆χ 2 = 1.78.
However, the freedom in σ ann v also means we can increase the Galactic centre Jfactor in order to lower the annihilation cross-section and still fit the excess. The tension between the dSph and the Galactic centre excess can therefore be relieved.

Conclusions
In this article we have presented global scans of combined dark matter-parton EFT operators for both real and complex scalar dark matter, including constraints from cosmology, indirect, and direct detection experiments. We have produced posterior distributions and profile likelihood maps of the model parameter spaces, and thus provided a state of the art, comprehensive, and -more importantly -coherent picture of scalar dark matter. We have investigate the prior dependency of the posteriors and seen that while the main features, i.e. peaks, are prior independent the tails of the distributions are decidedly not. We included a profile likelihood analysis for this very reason, and in doing so we have seen the complementary nature of marginalising and profiling. The marginalised posterior show us the bulk of probable parameter space while the profiled likelihood reveals small patches still allowed.
From these distributions we see that, of all current experimental results, the relic density has the greatest impact on the model parameters, as it is an actual measurement and not a limit. All the operator coefficients combined must be just right to provide the correct relic density. The most common configuration, and therefore the most probable, is when a single operator dominates and the others are weak. In contrast with the relic density the null results from LUX, the Fermi -LAT dSph analysis, and the Planck CMB anisotropy measurement have a somewhat limited impact on the operator coefficients. Instead, they do have effect on the dark matter mass, strongly disfavouring dark matter masses below 100 GeV.
The direct (LUX) and indirect (Fermi -LAT, Planck) searches each provide important complementary information. The direct detection experiment allows us to discriminate between operators that do and do not contribute to the spin-independent scattering crosssection. The indirect detection experiments similarily allows us to discriminate between p-wave and s-wave operators. In addition, the indirect detection measurements are the ones that constrain the m χ < 10 GeV region. A region which was allowed in earlier work (i.e. ref. [18]) that omits the indirect constraints.
Small patches survive in this low mass regime where the dark matter mass is around the charm or the bottom quark mass due to resonance-like effects in the relic density calculation. These regions could very well be excluded by the Fermi -LAT dSph data if the publicly available likelihood from Fermi -LAT would be extended to photon energies below 0.5 GeV.
We also considered the impact of including in our analysis the Galactic centre excess as signal of dark matter. Although standard astrophysical sources might explain the excess, such as unresolved millisecond pulsars, we showed the effect of including it in our likelihood to illustrate the effect of a measurement on the model parameters. The results are dramatic -a specific mass (40-60 GeV) and operator (R2/C2) are preferred to fit the excess. This region of the parameter space is however in tension with the absence of excess gamma rays from the Milky Way dwarf spheroidals in the real scalar case. So if one believe the Galactic centre excess is due to dark matter, then the quite general real scalar DM model would be excluded. The complex scalar candidate is able to avoid the tension by having a JHEP09(2016)077 lower annihilation cross-section and still have the correct relic density. This is achieved by combining a p-wave operator (C4) and a s-wave operator (C2 and/or C6.) Our toolchain is model-independent which means that, in principle, extending these types of combined EFT models -adding non-parton operators, or fermionic dark matteror looking a more complete models such as simplified models is as easy as writing down the Lagrangian. The main caveat being that adding more degrees of freedom without additional experimental constraints is inadvisable until more informative data become available. the "Fenomenologia y Cosmologia de la Fisica mas alla del Modelo Estandar e lmplicaciones Experimentales en la era del LHC" (FPA2010-17747) MEC projects. Finally also thanks to the Severo Ochoa MINECO project: SEV-2015-0398. This work was supported by Grant ST/N000838/1 from the Science and Technology Facilities Council (UK). R.T. was partially supported by an EPSRC "Pathways to Impact" grant. The work of TMPT is supported in part by NSF grant PHY-1316792 and by the University of California, Irvine through a Chancellor's Fellowship. We gratefully acknowledge the use of the LISA cluster (Amsterdam).
Open Access. This article is distributed under the terms of the Creative Commons Attribution License (CC-BY 4.0), which permits any use, distribution and reproduction in any medium, provided the original author(s) and source are credited.