Towards understanding thermal history of the Universe through direct and indirect detection of dark matter
Abstract
We examine the question to what extent prospective detection of dark matter by direct and indirect detection experiments could shed light on what fraction of dark matter was generated thermally via the freezeout process in the early Universe. By simulating putative signals that could be seen in the near future and using them to reconstruct WIMP dark matter properties, we show that, in a model independent approach this could only be achieved in a thin sliver of the parameter space. However, with additional theoretical input the hypothesis about the thermal freezeout as the dominant mechanism for generating dark matter can potentially be verified. We illustrate this with two examples: an effective field theory of dark matter with a vector messenger and a higgsino or wino dark matter within the MSSM.
a,b]Leszek Roszkowski, a,c]Sebastian Trojanowski d]and Krzysztof Turzyński Prepared for submission to JCAP
Towards understanding thermal history of the Universe through direct and indirect detection of dark matter

National Centre for Nuclear Research,
Hoża 69, 00681 Warsaw, Poland 
Department of Physics and Astronomy, University of Sheffield,
Sheffield S3 7RH, United Kingdom 
Department of Physics and Astronomy, University of California, Irvine,
California 92697, USA 
Institute of Theoretical Physics, Faculty of Physics, University of Warsaw,
ul. Pasteura 5, 02093, Warsaw, Poland
Keywords: dark matter theory, dark matter experiments, physics of the early universe
ArXiv ePrint: 1703.00841
Contents
1 Introduction
A detection of nongravitational interactions of dark matter (DM) particles would be one of the biggest breakthroughs in contemporary physics. There are many experimental strategies to search for such a signal, which include direct detection (DD) of DM particles by their scattering off heavy nuclei in deep underground detectors, indirect detection (ID) of products of DM annihilations or decays with space and ground based telescopes, and collider searches. So far, these efforts have led to stringent limits on DM interactions, but many of the favorable and very promising theoretical scenarios have not been fully probed yet, e.g., the prototypical scenario of a weakly interacting massive particle (WIMP), which we will denote by , such as the lightest neutralino in the Minimal Supersymmetric Standard Model (MSSM) (for review see, e.g., [1]).
Gravitational interactions of DM point towards a very precise value of its contribution to the energy density of the Universe, [2]. It is most often considered to have arisen in the socalled thermal production of DM in a freezeout process: the decoupling of the nonrelativistic DM particles that previously were in kinetic and chemical equilibrium with the thermal plasma, see, e.g., [3], and introduces a strict correlation between the mass and the interaction strength of DM particles. One may, however, also consider DM candidates that do not fit to the standard thermal production mechanism, but can still be effectively produced in other ways (for a review see, e.g., [4]). If the DM particles couple too strongly to the SM sector, the DM abundance is reduced and one needs to consider an additional mechanism of DM production.
In this paper, we aim at studying the properties of DM particles in a new way with the ultimate goal of understanding better the thermal history of the early Universe. We adopt a view that within the next decade or so collider searches for physics beyond the Standard Model will give negative results, but there will be a WIMP DM detection in one or more of the following experiments: Xenon1T [5] direct search for DM particles, the FermiLAT search for DM induced rays [6] and the Cherenkov Telescope Array (CTA) [7], a groundbased telescope which will also be looking for ray signal from DM annihilation and is expected to begin operating in some years from now. However, we are going to analyze the prospective signal cautiously and conservatively, taking into account known uncertainties. One can envisage a few possible scenarios. Given the sensitivities of the experiments, it is quite plausible that a positive DM signal would correspond to interactions strength larger than required for thermal DM production, which would point towards the existence of another DM component, presumably the particles of the same identity but originating from a process different from freezeout. Another goal of ours is to characterize, to the largest possible extent, that additional component and the mechanisms that gave rise to it. To this end, we simulate the signals that can be obtained by the collaborations in the case of a discovery and, subsequently, we use them to derive the basic properties of the DM particles, e.g., and the today’s value of the thermal average of the annihilation cross section times the Møller velocity (from now on referred to as annihilation cross section), , with their respective uncertainties, including astrophysical interpretation (for previous similar studies for the FermiLAT see [8, 9]). In our approach it is possible to discuss these constraints at various levels of theoretical assumptions: from a modelindependent level in which it is enough to introduce a set of phenomenological parameters to describe the relevant experimental signatures that we take into account, through an effective field theory of DM to very specific MSSM examples.
To this end, in Section 2 we perform a modelindependent analysis and conclude that, given the current limits on and known uncertainties, it would be very difficult, if not impossible, to rule out the possibility that DM consists only of the freezeout component, even in the simplest models, not to mention scenarios of multicomponent dark matter or with Sommerfeld enhancement (SE). In Section 3, we move on to describing DM within the framework of an effective field theory of a vector messenger, focusing on additional information that can be obtained from including positive signals from direct detection experiments. In Section 4, we analyze the specific case of the MSSM with neutralino DM. In this framework, we can discuss the mechanism for generating a nonthermal component of DM from decays of heavier particles such as a gravitino or an axino that belong to the category of extremely weakly interacting massive particles (EWIMPs). Since such particles were thermally generated at earlier stages of the evolution of the Universe, a determination of nonthermal components constrains the value of the reheating temperature. We conclude in Section 5. Technical details, in particular those regarding reconstruction of DM properties from the signals from direct and indirect detection experiments, are deferred to the Appendices.
2 General WIMP DM
We assume that WIMP DM could be produced in the early Universe both in freezeout, as well as in nonthermal processes, e.g., latetime decays of some heavier species. The WIMP DM relic density can then be schematically written
(2.1) 
where and are the contributions to the WIMP DM relic density from the freezeout and nonthermal processes, respectively.
The first term on the r.h.s. of Eq. (2.1) depends on the annihilation cross section of DM particles which can be often approximated as
(2.2) 
with coefficients and parameterizing the wave and wave contributions to , respectively. Eq. (2.2) can be used to relate the value of today which can be determined from the indirect detection of DM with the value of the annihilation cross section at the DM freezeout, , which in turn determines the value of . One can then make inferences about the limits on the additional DM component, presumably produced after thermal freezeout.
While the simple relation in Eq. (2.2) is not always valid, e.g., when coannihilations between the DM particles and other species are important around the DM freezeout [10, 11], one can still apply it for a conservative estimate, because additional interactions typically lead to an increase of the effective annihilation cross section at freezeout. There are, however, some scenarios in which these estimates do not apply. One possibility is that is smaller than due to interactions with other particles, e.g., resonant annihilation or Sommerfeld enhancement. One can also envision a more complicated DM sector, consisting of two or more stable particles that have undergone freezeout separately. Therefore, in the following we will first discuss simple scenarios in which (2.2) can be applied and then turn to a few representative models with the presentday annihilation crosssection larger than at freezeout. We also note that within a specific DM model, such as neutralino DM with a gravitino or axino companion in the MSSM investigated in Section 4, one does not need to resort to the approximation (2.2), as a precise relation between and is known.
2.1 The simplest case: singlecomponent DM with presentday annihilation crosssection not larger than at freezeout
The parameter space that we scan over for a modelindependent study is defined in terms of the WIMP DM mass , its ID cross section, , as well as the branching ratios that describe the final states for DM annihilation as shown in Table 1. We allow DM to annihilate into several final states or their mixtures. In particular we focus on , , and final states that correspond to distinct annihilation spectra. There are many other Standard Model final states, but usually they can be accurately approximated by one of the final states above or by the combinations thereof. However, we do not discuss here the purely leptonic final states, e.g., , since they are typically easily distinguishable from other cases that we take into account if is large enough. For a more detailed discussion about the reconstruction for different final states, see [12].
In order to remain as model independent as is possible, in this Section we do not assume any cross correlation between and . In practice, this is equivalent to excluding DD entirely from our considerations, because of the wellknown degeneracy in the recoil spectra for GeV. This is particularly relevant for our analysis, because it is just the mass range obtained by combining our assumption that is larger than the value required for thermal WIMP DM density with the current limits from null results of DM searches in dwarf spheroidal galaxies (dSphs) [13, 14] (see also [15] for the limits relevant for monochromaticlike spectra). In other words, irrespectively of whether DM DD signal is detected or not, for a given one can always adjust to account for that part of experimental input.
Symbol  Parameter  Scan range  Prior distribution 

WIMP mass  log  
Annihilation cross section  log  
(indirect detection)  
Branching ratios for various final states  modified Dirichlet distribution  
(branching ratios add to 1) 
Once we determine the CL intervals for the quantities describing the wouldbedetected WIMP DM particles, i.e., and , we can estimate the contribution to the DM density originating from freezeout, using Eq. (2.2). The most conservative estimate comes the assumption of wave dominance, , which translates to . Then the additional contribution to the relic density can be obtained from Eq. (2.1). Examples of results of such a reconstruction (red squares) are shown in the left panel of Fig. 1 for a TeV DM with the annihilation cross section equal to . We assume that the DM particles annihilate either purely into , , or into final state with admixture of . The mass reconstruction for the pure hadronic final state is quite poor, because the corresponding annihilation spectrum is often degenerate with the spectrum obtained for other pure or mixed final states for different DM masses. However, this can be drastically improved if either small admixture of the monochromatic line is present in the spectrum (from or final state; magenta dots) or the final state is not purely hadronic but also contains a leptonic contribution as can be seen for the case with final state (blue triangles). On the other hand, it is important to note that the quality of the reconstruction of is rather sensitive to the derived values of and therefore even a poor reconstruction of the mass parameter can lead to constraints on the nonthermal contribution to the DM relic density.
However, if the DM particles interact more weakly, the quality of reconstruction decreases drastically leading to practically no conclusive results about the freezeout origin of DM for (green squares) and (grey squares) at a modelindependent level. We will later see that this can be improved, e.g., if a correlation between and is derived for a specific model of the DM interactions.
Nonzero values of translate to , which suppresses the abundance of WIMP DM at freezeout and therefore enhances . This is illustrated in the right panel of Fig. 1 for two values of and that lead to a nonnegligible or even dominant wave contribution to the annihilation cross section in Eq. (2.2) around the freezeout temperature .
These limits can be presented in a way that can be useful for future studies on the additional, nonthermal, DM production in the context of specific models. In Fig. 2, we identify the values of and of some underlying model (benchmark point) for which one can infer that a nonthermal component to the WIMP DM density is necessary. The region marked ‘unclear’ corresponds to the benchmark points for which the 95% CL range of the reconstructed values contains the value of the annihilation cross section consistent with the total WIMP DM density originating from freezeout. We note that in a purely modelindependent approach it is practically impossible to confirm the presence of nonthermal DM component, as large values of that would allow ruling out the possibility of freezeoutonly DM are in conflict with ID searches except for tiny corners of the parameter space. However, it is important to remember that the published limits on , e.g., those by the H.E.S.S. collaboration [16], depend on the assumed DM profile which we allow to vary in our reconstruction (see Appendix A.3).
With a specific choice of the annihilation final state, we can discuss the results in a finer way. In Fig. 3 we show the minimum value of the reconstructed in the 95%CL region as a function of the benchmark point WIMP DM mass, , for several choices of the annihilation cross section for the benchmark point, . The limits are derived assuming that DM annihilates either purely into one of the final states that we consider or there is small admixture of the branching ratio into final state that leads to a monochromatic contribution to the spectrum. In the simplest case in which the DM annihilation in the early Universe is dominated by the wave contribution and coannihilations play negligible role, the minimum value obtained in the reconstruction can be directly translated into a lower bound on and therefore into limits on the contributions to the DM relic density can be derived. This is shown on the scales on the righthand sides of the plots in Fig. 3, in which the minimal value of required to obtain the correct DM relic density is shown. These values can also be viewed as conservative limits on since both nonnegligible wave contribution, as well as additional coannihilations typically tend to increase above the current value . The annihilation cross section needs to be at least about an order of magnitude above the freezeout value in order to allow the reconstruction of , while for lower values one would not be able to determine if the additional contribution to the DM relic density is needed. However, as we will see in the next section, this can be improved within a framework of a specific model.
As can be seen in the left panel of Fig. 3, the quality of reconstruction of for final state can be improved if a small admixture of a monochromatic line is present in the DM annihilation spectrum. Such an improvement is less pronounced for both and final states shown in the left panel of Fig. 3. In particular, in the case of this is so because of the characteristic monochromaticlike spectral feature that is already present in the annihilation spectra. It comes from splitting with the emission of soft and greatly increases the accuracy of DM mass determination. For this reason we do not show the lines with small admixture in the right panel of Fig. 3. In the plots we also show FermiLAT and MAGIC exclusion lines based on [13]. It is important to note that these lines in Fig. 3 cannot be directly compared to the usual exclusion lines shown in [13], because in our case they correspond to the benchmark values of the annihilation cross section, , which are larger than the minimum values of the reconstructed cross section, . It is, however, the latter quantity that is shown on the vertical axes of both plots in Fig. 3, so the limits that we show in Fig. 3 should be then rather thought as the lower limits on the DM mass for a given value of indicated in the plots.
The results presented above can be also used in the presence of invisible annihilation channels. Effectively, this is equivalent to having smaller which leads to a poorer reconstruction of DM properties.
2.2 More complicated examples
2.2.1 Multicomponent thermal DM
Scenarios of multicomponent thermal DM pose an additional difficulty for extracting DM properties from the DM ID data even if (2.2) is satisfied for each component separately. This is due to the fact that the DM ID rate is proportional to , where runs over all DM components and is the local density of the th component, which in turn is proportional to the averaged annihilation crosssection for this component at freezeout, see Eqs A.1 and A.2. It is therefore possible that these crosssections are larger than in the singlecomponent (purely) freezeout case and can be, in principle, detectable in the forthcoming DM ID experiments.
The prospects of determining DM properties are, however, significantly worse for these models compared to the case of singlecomponent DM with both freezeout and nonthermal contribution, because the DM ID rates scale as , where is the presentday local DM density. We illustrate this in the left panel of Fig. 4, where we show the reconstructed value of the annihilation cross section multiplied by the squared ratio of the local DM density of each DM component, , with respect to the total local DM density, . The benchmark point is the , scenario discussed in Section 2.1, with . Light blue squares correspond to a scenario in which only one DM component leads to a ray signal in DM ID. We assume that this component has the same mass and crosssection as the benchmark model, but now only a thermal relic density, i.e., . The green squares represent the reconstruction of the model in which there is one more DM component with thermal relic density and , while (wave) is chosen such that the total DM relic density of both thermal DM components saturates the Planck data. In this case the observed ray signal is a joint signal from both these DM particles. In both examples, the quality of reconstruction is much worse than for the benchmark model alone; it resembles that of the singlecomponent DM with a much smaller annihilation crosssection.
2.2.2 Sommerfeld enhancement
A significant modification of the DM annihilation crosssection occurs in the nonrelativistic regime when DM couples to a light mediator. This effect, known as Sommerfeld enhancement [17, 18, 19, 20, 21], invalidates the simple relation (2.2) and we will discuss it separately. Here, we approximate the Yukawa potential by the Hulthén potential which allows to calculate the enhancement factor analytically [22, 23] (see also [24, 25] for further discussion). We follow [26] to regularize velocity near a zero energy bound state and to describe the mediator sector with two parameters: the coupling constant between the DM particle and light mediator, , and the mediator mass, . When performing reconstruction, we allow them to vary in ranges and , respectively.
In the presence of the SE, the present value of the annihilation cross section, , can be larger than , which means that the nonthermal contribution to DM density is smaller or not necessary at all. This is illustrated in the right panel of Fig. 4, where the effects of SE for the benchmark model of Section 2.1 are shown: we plot the reconstructed value of the annihilation cross section around freezeout for a fixed value of presentday annihilation crosssection. Light green squares correspond to the assumption that the observed DM density has a fully thermal origin, while for light blue squares this requirement has not been imposed. We can see that the effect can be large or small, depending on the SE factor determined by and ; in particular we can obtain the same BP for DM ID as discussed above, but with the around freezeout at the thermal level for and . We conclude that in the presence of SE it is impossible to infer whether a prospective observation of DM ID in forthcoming experiments signifies the need for nonthermal contribution to DM density. This is hardly surprising, as the MSSM wino is commonly known to be ruled out as a thermal DM candidate event with the current DM ID sensitivity [28, 27, 29].
Once the additional constraint on the thermal relic density is imposed, one expects the reconstruction of the annihilation cross section around freezeout to be much better, as illustrated by light green points in the right panel of Fig. 4. However, due to the freedom in obtaining large SE, the interplay between this strong constraint and the reconstruction based on DM ID signal is only via the reconstruction of the DM mass. Note that would remain practically unconstrained if only the DM relic density constraint was considered without DM ID signal taken into account.
3 The effective field theory of DM: an example
The quality of reconstruction of the annihilation cross section can be improved if, in addition to the signal coming from ID, DM is discovered in DD experiments thanks to possible correlations between scattering amplitudes in both types of searches [30]. Such correlations are, however, never ideal as typically is governed by DM couplings to heavy fermions, e.g., or , or gauge bosons, while DD rates are associated with DM couplings to light and quarks. Although limited, such correlations may improve the reconstruction within the framework of a given model, because fitting the same ID signal with different final states would require a shift in and, correspondingly, in , and the latter might be in tension with the simultaneous fitting of the DD signal. Such improved reconstruction of the annihilation final state can lead to better limits on the annihilation cross section and, subsequently, on the nonthermal contribution to the DM relic density.
We illustrate this effect within a framework of Effective Field Theory (EFT) approach to study DM couplings [31, 32, 33, 34, 35, 36]. In particular, we assume vectorlike couplings between DM and the SM particles (see, e.g., [37, 38, 39, 40, 41, 42, 43]). As we focus on the heavy DM particles, with TeV, even with prospective ID and DD signals we typically remain beyond the reach of the LHC [44, 45] (for EFT approach to DM DD see, e.g., [46, 47]).
We assume that the DM particles have vectorlike effective couplings to the SM chiral eigenstates, , where corresponds to the cutoff scale of the UV completion of the model. The hierarchy between the coefficients (also assumed), is such that DM couples dominantly to 3rd generation fermions. In the mass eigenstate basis this leads to the following Lagrangian that contains both vector and axialvector couplings
(3.1)  
The Lagrangian in Eq. (3.1) can be written in the chiral basis
(3.2) 
with coefficients
(3.3) 
In this effective model DM can then annihilate into one of the following final states or their mixture: , , and . The annihilation cross section would have both the vector and axialvector contributions [49, 50]. Typically wave contributions to are suppressed, but, for completeness, we take them into account (in addition to the dominant wave contributions) both when treating DM ID, as well as for the calculation of the DM relic density.
Although a direct coupling of the DM particles to and quarks is absent in Eq. (3.1), the respective Wilson coefficients will be generated once the Renormalization Group (RG) evolution is taken into account from the renormalization scale down to the nuclear scale GeV [48, 51]. These Wilson coefficients for light and quarks determine the DD scattering cross section, . The absence of direct couplings of DM to and quarks at the tree level allows us to obtain both and that are large enough to be detected for the purpose of our reconstruction, but at the same time not yet excluded.
Parameters  Ranges 

TeV  
TeV  
, , , , 
We vary coefficients , , , and , as well as the scale and the DM mass, , as shown in Table 2. The reconstruction is performed for benchmark points characterized by the dominant annihilation channels into or . The DM mass, TeV, and the annihilation cross section, , are chosen at the same level as in the modelindependent analysis. The () benchmark point corresponding to the largest value of is characterized by (). Our typical values of for benchmark points, as well as in the CL reconstructed regions, exceed by a factor of a few. The quality of the EFT approach in such a case has been studied, e.g., in [40] for axialvector coupling: in such a scenario the EFT approach is either valid within some error or may lead to underestimated values of cross sections. In the latter case our effective model would correspond to the conservative lower limit on nonthermal contribution to the DM relic density which could be further improved if mediator is included in calculations.
In the left panel of Fig. 5 we show the correlation between and for points in the considered EFT model with TeV and varying annihilation branching ratio to final state. In particular, the points with (green circles) shown strong correlation. This is not surprising, since for the fixed mass and annihilation cross section, both and effectively depend on a single quantity, i.e., the cutoff scale . On the other hand, once one allows more efficient annihilation into different final states, the correlation becomes weaker.
We can now see how an interplay between DD and ID experiments can help to improve reconstruction of the DM properties. In the large mass regime, the DD experiments have no sensitivity to DM mass, as discussed in Appendix A.2, but if the DM signal is observed in the ID experiment one can obtain relatively accurate values of and the dominant annihilation branching ratio. This can lead to a stronger correlation in expected signal between ID and DD experiments, which, in turn, results in an improved reconstruction of both and .
Finally, one can translate the improved limit on into the limit on the nonthermal contribution to the DM relic density. We show such an improved reconstruction of in the right panel of Fig. 5 (compare with the left panel of Fig. 1 for which we used the same and for the benchmark points). Compared to the analysis presented in Section 2, we note an improved mass reconstruction in the ID experiments with respect to the similar analysis discussed in the modelindependent section, which comes partly from a smaller set of final states in the considered EFT model (specifically, the lack of and final states). On the other hand the improvement in the reconstruction of is to a large extent the effect of the correlation between the ID and DD experiments.
4 Neutralino DM
In the preceding Sections, we have shown what we can learn about the generation of WIMP DM without invoking specific assumptions about the underlying microscopic model of DM. We then discussed how such a reconstruction can be improved by adding simple assumptions about the correlations between the DD and ID rates. We illustrated the latter feature within a framework of a simple EFT model. In both cases we focused solely on DM itself, i.e., we ignored any possible additional contribution to from coannihilations that could become important if the DM particles are massdegenerate with some other species in the thermal plasma in the early Universe.
By taking into account coannihlations one allows in general at least one more degree of freedom associated with the mass difference between and the heavier species. This, in principle, could make it much more difficult to draw any conclusions about the contributions to the DM relic density from reconstruction of the current value of the annihilation cross section, , unless some other nontrivial relations between all these quantities are present in the model. We will illustrate this point by performing a study of a prototypical example of a WIMP DM particle, i.e., the lightest neutralino, , which appears in the context of supersymmetric theories.
4.1 Reconstruction of nonthermal component of DM with wino and higgsino DM
In general, the lightest neutralino is a linear combination of the superpartners of the gauge bosons and the Higgs bosons. In our study we focus on scenarios in which the neutralino is higgsino or wino dominated. It is known that for these type of DM particles freezeout can provide a correct amount of DM density for about TeV and TeV, respectively, while for lower DM mass one obtains too low , so a discovery of such a WIMP DM particle at a lower mass scale would signify the necessity of the nonthermal contribution to (or another DM species).
We also note that in both scenarios coannihilations play an important role in establishing the DM relic density in the early Universe. In particular, the difference in mass between and the lighter chargino, as well as second to the lightest neutralino in the case of higgsino DM, is typically bound to no more than several GeV. For such a small mass difference, the DM freezeout density is not very sensitive to the precise determination of the mass degeneracy, but it depends on the model parameters that also determine the composition of the lightest neutralino which in turns settles the actual values of both and . Such a set of nontrivial correlations between DD and ID rates, as well as the DM relic density allows for successful reconstruction of the nonthermal component of the DM relic density even in the presence of coannihilations during thermal freezeout.
We present our results for neutralino DM for two benchmark points that correspond to higgsino DM with GeV and (predominantly) wino DM with TeV. The cross sections for direct and indirect detection for both benchmark points have been chosen to lie just below the current upper bounds obtained by the LUX Collaboration [52], FermiLAT + MAGIC [13, 14] and H.E.S.S[16], respectively. It is important to note that these bounds should not be treated as a strict limits since they correspond to the standard choice of the astrophysical parameters while we allow them to vary in the scan (see Table 3).
We perform a numerical analysis of a 10parameter phenomenological MSSM; details of the numerical procedure and the applied constraints are described in the Appendix B. The CL regions obtained after fitting to the particle physics constraints and to the DM signal mock data set are shown in Fig. 6 for both of our benchmark points. We present the results in the plane where . The quality of reconstruction allows to constrain the necessary nonthermal contribution to the neutralino DM relic density, which varies between and for the higgsino DM case, and between and for the wino DM scenario. The upper limits of the reconstructed CL regions correspond to smaller freezeout contribution, which can be roughly translated to larger , implying that the current constraints from DM ID become important. The lower limits on are driven by the quality of reconstruction of the DM ID signal, which is mainly affected by astrophysical uncertainties: a change in the nuisance parameters described in Table 3 can make a model with smaller mimic the DM ID signal. For both our benchmark points, the CL regions would be significantly smaller if one neglected the astrophysical uncertainties.
The absolute quality of the DM mass reconstruction is similar for both benchmark points, irrespective of . This follows from a more pronounced monochromaticlike spectral feature coming from final state for . On the other hand, is better reconstructed for the higgsino DM benchmark point with GeV. It is mainly because of the interplay between the DM DD and ID which limits the impact of varying . As can be seen from Eqs (A.1), (A.2) and (A.4), modifying in such a way that remains unchanged requires a change in such that is roughly constant. For a given this can be achieved by adjusting the gaugino component of the higgsino. In principle, there are two adjustable gaugino mass parameters, i.e., the bino mass, , and the wino mass, , which appears to give enough freedom to simultaneously fit ID and DD signals for modified . However, increasing the bino component, would lead to a decrease of the annihilation cross section and an increase of . Therefore, varying the bino component does not help in ‘hiding’ the DM signal in both DD and ID; since for a varying wino component the DD and ID constraints are nondegenerate, the DM reconstruction becomes more accurate in the higgsino DM case.
4.2 Origin on nonthermal DM and constraints on the reheating temperature
In the previous section we studied the quality of reconstruction of the nonthermal contribution to the DM relic density of higgsino and wino DM without specifying its actual origin. Interestingly, within the MSSM there exists a mechanism for the production of such additional DM particles; it involves extremely weakly interacting massive particles (EWIMPs), with interactions so weak that EWIMPs decay to DM particles long after the freezeout. Theoretically motivated EWIMPs include the axino [53, 54, 55, 56, 57] (see also [58] and for a review [4, 59]) and the gravitino [60, 61, 62, 63, 64, 65] that appear in supersymmetric models. Both EWIMPs have very low interaction rates, suppressed by a high cutoff scale, but nonetheless they can be effectively produced in, e.g., scatterings and/or decays of other particles that were in thermal equilibrium themselves [57, 65, 66, 67, 68, 69, 70, 71].
EWIMPs interact too weakly to be detected in any of the current or foreseeable future experiments searching for particles beyond the SM. Typically, the limits on these particles come from cosmological considerations about the Big Bang Nucleosynthesis (BBN) [72, 73, 74, 75], Large Scale Structure (LSS) formation [76] and Cosmic Microwave background (CMB) radiation [77, 78]. However, within a given model, one can look for collider signals associated with species heavier than EWIMP, see, e.g., in [79, 80, 81, 82, 83, 84, 85, 86, 87, 88] (for more recent studies related to this topic see, e.g., [91, 92, 93, 94, 89, 90]). Alternatively, additional information can be obtained from DM ID signals in the scenarios of decaying gravitino or axino DM [98, 95, 96, 97, 99, 100]. Altogether, these constraints on DM particles and their EWIMP companions translate to constraints on the thermal history of the Universe.
We will illustrate this issue for the neutralinos that come from latetime decays of heavier axinos or gravitinos after neutralino freezeout. In that case, eq. (2.1) can be written as
(4.1) 
where denotes the yield of EWIMPs from socalled thermal production of EWIMPs,^{1}^{1}1Thermal production of EWIMPs should not be confused with freezeout thermal production of WIMPs. In the latter case one assumes that WIMPs are in thermal equilibrium in the early Universe. This is, however, not true for EWIMPs for which one only assumes that they originate from decays and scatterings of other species that are in thermal equilibrium., , where is the presentday entropy density and is the critical density, and we assume . For the gravitino or the axino, their production from thermal plasma is governed by nonrenormalizable operators and therefore the EWIMP yield is sensitive to the reheating temperature, , which can be identified with the maximal temperature in the radiation dominated (RD) epoch in the evolution of the Universe.^{2}^{2}2For a recent discussion of corrections to the axino yield due to noninstantaneous reheating period see [101]. In the case of the renormalizable operators the final yield of EWIMPs would be proportional to the annihilation cross section [102], , as in any generic freezein mechanism [103].
4.2.1 Axino EWIMP
Axinos can be produced effectively if the reheating temperature of the Universe after inflation, , is smaller than the temperature at which these EWIMPs establish thermal equilibrium or they could otherwise easily overclose the Universe. The amount of produced axino EWIMPs is typically proportional to , unless the reheating temperature becomes small and the Boltzmann suppression of the number densities of the involved particles needs to be taken into account. For values of smaller than the freezeout temperature one also needs to consider a modified evolution of the Universe around freezeout. We take all this effects into account following [101, 104]. We consider the supersymmetrized version of KimShifmanVainsteinZakharov type of axion models [105, 106]. We also assume that the axino is lighter than the gluino, . This is merely a technical assumption: otherwise, depending whether the gluinos decay into neutralino DM before or after freezeout, this extra contribution to DM is either erased or stays intact, respectively [107]. As a result, we obtain one of two simple cases, with the choice sensitive to the details of the MSSM spectrum. With this assumption, the axino mass can be varied between and without affecting much the nonthermal contribution to the neutralino relic density which depends essentially only on for a given point in the parameter space of the MSSM. Therefore, we can present previously derived bounds on as constraints on as shown in Fig. 7. As expected the reheating temperature in this scenario is confined to low values and the order of magnitude of its value can be determined for both our benchmark points. The uncertainty on is mainly driven by the lack of accuracy in determination of and , but to some extent it also comes from the fact that the production of axinos for such low is sensitive to the rest of the supersymmetric spectrum.
4.2.2 Gravitino EWIMP
Gravitino interactions are suppressed by the Planck mass, GeV, which is several orders of magnitude larger than the PecceiQuinn scale relevant for the axino, GeV [108]. As a result in order to obtain a correct amount of nonthermally produced neutralinos coming from latetime decays of gravitinos one needs to consider significantly larger reheating temperatures. Importantly, the suppression by also results in longer lifetime of the gravitino that could lead to latetime decays violating the predictions of the BBN [109, 110]. One also needs to take into account possible implications of the fact that the neutralinos produced in gravitino decays may form a warm fraction of neutralino DM. In particular, it has been shown that this effect can become important for wino DM with mass around GeV [111].
Here we assume that the gravitino is significantly heavier than the neutralino, . More precisely, we assume that TeV [112], which means that the gravitino is also significantly heavier than all the gauginos in our framework. The yield of thermally produced gravitinos is then practically independent of the gravitino mass. Similarly to the axino, the gravitino cannot be arbitrarily heavy in order to avoid a scenario in which it decays into neutralino before neutralino freezes out and there is no nonthermal contribution to the neutralino DM density.
The nonthermal contribution to the neutralino DM density in our scenario reads [65, 67, 68]
(4.2) 
Unlike in the axino case, for which there was no such a simple estimate of like the one given in Eq. (4.2), because for low values the axino yield is not simply proportional to the reheating temperature and depends quite significantly on the MSSM spectrum [101], it is now straightforward to employ Eq. (4.2) to translate the limits on into constraints on . We illustrate this in Fig. 8 for both our benchmark scenarios. As expected the reheating temperature required to obtain the correct relic density is now significantly larger than for the axino. The quality of reconstruction is limited by both the uncertainty of and , but in both scenarios can be determined with very good accuracy for wino DM benchmark scenario and up to a factor of for the higgsino DM case due to larger uncertainty on the DM mass associated with the characteristic monochromaticlike feature present in annihilation spectrum for final state that is more pronounced for larger DM masses.
5 Conclusions
In this paper, we addressed the question whether prospective detection of DM in the next generation of experiments could shed light on whether DM was generated thermally in the freezeout process in the early Universe. To this end, we simulated signals that could be seen in the indirect detection experiments FermiLAT and CTA, as well as in the direct detection search in the Xenon1T experiment for DM particles with annihilation cross section significantly larger than required for thermal freezeout. We then reconstructed the mass and the annihilation cross sections of such DM particles from these wouldbe signals and checked if the nonthermal component could be further understood in terms of appropriate underlying models.
We showed that in the modelindependent approach the answer is negative except for a thin sliver in the parameter space assuming an wave annihilation. It could even become not possible once one considers, e.g., the Sommerfeld enhancement of the annihilation cross section or multicomponent DM scenarios. Hence some theoretical input is required for distinguishing the additional nonthermal component in the next decade or so. We discuss two such theoretical scenarios, varying in the degree of complexity: an EFT of DM with (axial) vector messenger and the MSSM. We showed that even with rather general assumptions about the EFT, the reconstruction of DM properties improves to the extent that the nonthermal component can be identified from the reconstruction of prospective ID and DD signals. We then turned to the MSSM and showed examples of benchmark points for which and/or can be reconstructed even more precisely – either because DM particle annihilate to with important contribution from the Sommerfeld enhancement or because DM DD is sensitive to the gaugino composition of the neutralino. In addition, nontrivial correlations between the DD and ID rates, as well as the DM relic density allows to infer conclusions about the nonthermal contribution to the DM relic density in such a scenario even in the presence of effective coannihilations around the DM freezeout.
This enhanced precision of the determination of the nonthermal component can be used to study processes from which this component originated. We demonstrated this with two examples of EWIMPs constituting this component: the axino and the gravitino. In each case we determined the corresponding reheating temperature and concluded that – compared to the present particle physics constraints – future prospective detection of DM signal would lead to a moderate improvement of the allowed range of this parameter for the axino EWIMP, while for the gravitino EWIMP it would significantly narrow down the allowed range.
Taken together our results illustrate the possibility that the discovery of the DM signal alone – though a qualitative step forward in characterization of the composition of the Universe – is likely to give rise to a number of further questions. Answering these questions could require working in a specific theoretical framework, which would have to be inferred from other experiments.
Acknowledgments
ST would like to thank A. Hryczuk, T.M.P. Tait and P. Tanedo for helpful comments. LR is supported in part by the National Science Council (NCN) research grant No. 201518AST200748 and by the LancasterManchesterSheffield Consortium for Fundamental Physics under STFC Grant No. ST/L000520/1. ST is supported in part by the National Science Centre, Poland, under research grant DEC2014/13/N/ST2/02555, by the Polish Ministry of Science and Higher Education under research grant 1309/MOB/IV/2015/0 and by NSF Grant No. PHY1620638. KT is partly supported by a grant 2014/14/E/ST9/00152 (National Science Centre, Poland). The use of the CIŚ computer cluster at the National Centre for Nuclear Research is gratefully acknowledged.
Appendix A Reconstruction of general DM particles
a.1 Gamma rays from DM annihilations
Most of the rays that reach detectors have an astrophysical origin which needs to be carefully taken into account when looking for signal from DM annihilation. However, one expects that DMinduced rays are typically produced in dense regions with large and therefore they should manifest themselves as a small directional excess over the isotropic background. One such obvious region of interest for DM searches is the Galactic center (GC) where a significant increase of ray flux from DM annihilations is expected due to the peak in the DM mass distribution. This is in particular true for the original NFW profile [113], as well as for the Einasto profile [114, 115], while for the DM distributions with cores in the GC, e.g., the Burkert profile [116], the increase of DMinduced ray flux is not so pronounced.
Other promising sources of the rays from DM annihilations are dSphs around the Milky Way. They are some of the most DMdominated nearby objects known and therefore the DM searches focusing on dSphs suffer less from astrophysical background than in the case of the GC. However, the corresponding factors are typically smaller and their exact values are more difficult to determine. This introduce additional challenge when interpreting the experimental results, but still the most stringent limits on up to date come from null DM searches in the combined analysis of FermiLAT and MAGIC observations of dSphs [13, 14].
The total differential flux of the rays from a source with the angular size is given by
(A.1) 
where stands for the ray spectrum for the annihilation final state with the corresponding branching ratio , the total annihilation cross section is given by , while the second term in the sum corresponds to the contribution from the Inverse Compton scatterings of the DMinduced electrons from the GC (see, e.g., [117] for more details). In the case of the dSphs one does not expect significant contribution from the secondary ray emission due to lower energy densities of photons than in the GC [118].
The astrophysical uncertainties connected to the DM distribution, , are hidden in the factor which is defined as the integral of along the line of sight (l.o.s.) for the corresponding angular size of the DM source
(A.2) 
In the case of the GC, we employ the generalized NavarroFrenkWhite (NFW) DM profile for that reads [119]
(A.3) 
where we fix and the distance of the Solar System from the GC, , while we let the remaining parameters, i.e., and , to vary freely in the scan. For the dSphs we use the factors and their uncertainties from [120].
a.2 Direct detection of DM particles
Another techniques of WIMP DM searches employs its possible scattering off heavy nuclei in deep underground detectors. The signal induced by DM particles can then be recorded as light coming from scintillation photons produced by deexciting atoms, charge connected to atomic ionization or heat from phonons in crystals. In order to improve discrimination of the signal and background events coming from electron recoils one often employs two of the aforementioned detection strategies in a single experiment (for review see, e.g., [121]).
The distribution of the number of expected events, , as a function of the recoil energy, , of the DM particles off nuclei is given by
(A.4) 
where is the spinindependent WIMPproton scattering cross section, is the atomic number of the nuclei, the reduced mass of the system is given by , is the nuclear form factor, while the integral goes over the relative velocity between the WIMP DM particle and the nucleus, , and the distribution of the WIMP velocities, , has a cutoff at the galaxy escape velocity, . The minimum velocity that can result in an event with the recoil energy is given by , where is the mass of the nucleus and is the reduced mass of the nucleus system. We assume that the DM particles scatter elastically off nuclei and that there are no important isospin violation effects. Inelastic scatterings are typically suppressed with respect to the elastic ones [122, 123] and we neglect such possibility in our study (for reconstruction with inelastic scatterings and/or isospin violation effects taken into account see, e.g., [124, 125, 126]).
If the WIMP DM mass is larger than the nucleus mass, the minimal velocity, , is insensitive to and therefore the recoil energy distribution is simply proportional to as shown in Eq. (A.4). This makes it impossible to reconstruct the properties of such heavy DM particles from only DD experiments at a level of modelindependent WIMP DM since only the ratio can be inferred from the data. However, the interplay between DD and ID can help to overcome this problem [127, 128, 129, 12], especially in a framework of particular models in which one can benefit from the cross correlation between and .
a.3 General methodology of reconstruction
In order to determine the quality of reconstruction of the DM properties, we first generate a signal mock data set for wouldbe discovered DM particles defined by our assumed benchmark points. The reconstruction is then performed by scanning over the parameters of the models that we consider (see Tables 1, 2 and 4), as well as over the nuisance parameters that define the astrophysical uncertainties (see Table 3). For each point in the scan we estimate the quality of fit of the obtained data set to the fixed signal mock data set of the benchmark point by evaluating the respective likelihood functions. We finally present our results in 2dimensional plots showing confidence level (CL) regions defined by the condition .
Symbol  Parameter  Scan range  Prior distribution 

Circular velocity  km/s  Gaussian  
Escape velocity  km/s  Gaussian  
Local DM density  Gaussian  
NFW slope parameter  Gaussian 
The details of the likelihood functions and the reconstruction methodology can be found in [12]. In our reconstruction we base on the code developed for that study. Here we only briefly summarize the main aspects of the procedure. In the case of the GC we analyze ray signals from DM annihilations that will be observed by the CTA for which we assume hours of observational time. We take into account the background from the Galactic Diffuse Emission (GDE) and cosmic rays (CRs). In the fitting, we use the binned Poisson likelihood function convoluted with the Gaussian functions describing the uncertainties of the background estimation, where the bins are defined both in the photon energy, as well as for the four spatial regions in the sky around the GC. In the study of the DMinduced ray signal from the dSphs we focus on DM searches performed by the FermiLAT. For this analysis, we assume years of observation and an extended set of dSphs. The likelihood function in the fitting procedure follows [120] from where we also take the uncertainties on the factors.
In the case of the DM DD we employ a binned Poisson likelihood in which the expected signal is obtained by the integration of Eq. (A.4) in each energy bin. The residual background in DD experiments is typically small due to a high purity of target material and good separation from external sources of radioactivity. In addition, one of the main sources of the background, namely the electronic recoils, can be discriminated from the nuclear recoils thanks to a combination of two detection strategies, as in the dualphase Xenon1T time projection chamber that records signal both as scintillation light and ionization charge.
Appendix B Reconstruction of neutralino DM particles
b.1 Direct detection and gamma rays for indirect detection of neutralino DM
Higgsino and wino DM particles annihilate dominantly into , , and final states with possible admixtures of other channels that involve nonSM particles, e.g., and . This leads to a continuous spectrum of photons originating from subsequently produced cascades of particles. Additional contributions to ray spectrum comes from internal bremsstrahlung processes, as well as, more importantly, from loopsuppressed, but nonnegligible annihilations into and pairs that result in monochromatic lines. For completeness, we also take into account secondary emission of rays from the inverse Compton scattering of DMinduced electrons, although it typically plays a subdominant role for higgsino and wino DM due to suppressed annihilation branching ratios into leptons. When generating the annihilation spectra we employ tables provided in [117, 130], with the additional contribution from secondary rays that follows [131]. To obtain the cross sections for monochromatic lines, as well as ray spectrum from internal bremmstrahlung processes we employ MicrOMEGAs 4.3.1 [132]. In the case of the annihilation channels with products from beyond the SM, e.g., or , which are not treated in [117], we employ SUSYHIT [133] to determine their branching ratios into the SM particles. We then use for these SM decay products the ray spectra from [117] obtained for the appropriately shifted CM energies. Such a simplified approach is the most relevant for the annihilation products produced at rest, which we find applicable to the points in our scan that have significant branching ratios into the nonSM particles. In other scenarios the approach is based on the assumption that the a priori smooth distribution of energies of final SM particles in the CM frame of the annihilating DM particles can be approximated by the average values of these energies. We would like to stress here that this approximation does not affect much our final results since most of the points have dominant annihilation final states into the SM particles. In the case of wino DM, one also needs to take into account the Sommerfeld enhancement which can significantly modify the annihilation cross section both when calculating the DM freezeout density and [19, 134] (see also [135, 27, 28, 29, 136, 137] for more recent studies). Here we follow [138, 139].
The dominant contribution to the spinindependent cross section, , for higgsino and wino DM comes from a channel Higgs boson(s) exchange unless is a pure state. In the latter case both the channel squark exchange, whose rate is suppressed by the squark mass, as well as loopinduced processes become more important [140], but the scattering cross section that one obtains in such a scenario typically lies below the experimentally accessible limits. In accord with our approach of focusing on the scenarios that can be probed by DM DD in the following years, we choose our benchmark points to have subdominant but nonnegligible mixing between gauginos and higgsino.
b.2 Methodology of reconstruction for neutralino DM
We generate the signal mock data set for both benchmark points and subsequently attempt to fit it with the signal corresponding to other neutralino DM scenarios obtained in the scan over the parameter space of the MSSM. We perform the scan over a 10parameter version of the MSSM with the parameters and their ranges shown in Table 4. We employ Mutlinest [141] for sampling the parameter space of the model. We use SOFTSUSY3.4.0 [142] to generate the mass spectrum and we take into account the requirement that the mass of the lightest SMlike Higgs boson, , fits the observed value [143] with GeV theoretical error. We also take into account the LHC limits on gluino and squark masses following [144], as well as basic constraints related to physics, for which we use SuperIso v3.3 [145], although they play minor role in our discussion. The constraints imposed in the scan other than fitting to the signal mock data set are summarized in Table 5.
Parameter  Range 

bino mass  
wino mass  
gluino mass  
stop trilinear coupling  
stau trilinear coupling  
sbottom trilinear coupling  
pseudoscalar mass  
parameter  
3rd gen. soft squark mass  
3rd gen. soft slepton mass  
1st/2nd gen. soft squark mass  TeV 
1st/2nd gen. soft slepton mass  GeV 
ratio of Higgs doublet VEVs  
Nuisance parameter  Central value, error 
Bottom mass (GeV)  (4.18, 0.03) [146] 
Top pole mass (GeV)  (173.21, 0.87) [146] 
References
 [1] S. P. Martin, A Supersymmetry primer, Adv. Ser. Direct. High Energy Phys. 21 (2010) 1 [Adv. Ser. Direct. High Energy Phys. 18 (1998) 1] [hepph/9709356].
 [2] P. A. R. Ade et al. [Planck Collaboration], Planck 2015 results. XIII. Cosmological parameters, Astron. Astrophys. 594 (2016) A13 [arXiv:1502.01589 [astroph.CO]].
 [3] E. W. Kolb and M. S. Turner, The Early Universe, Front. Phys. 69 (1990) 1.
 [4] H. Baer, K. Y. Choi, J. E. Kim and L. Roszkowski, Dark matter production in the early Universe: beyond the thermal WIMP paradigm, Phys. Rept. 555 (2015) 1 [arXiv:1407.0017 [hepph]].
 [5] E. Aprile et al. [XENON Collaboration], Physics reach of the XENON1T dark matter experiment, JCAP 1604 (2016) no.04, 027 doi:10.1088/14757516/2016/04/027 [arXiv:1512.07501 [physics.insdet]].
 [6] W. B. Atwood et al. [FermiLAT Collaboration], The Large Area Telescope on the Fermi Gammaray Space Telescope Mission, Astrophys. J. 697 (2009) 1071 [arXiv:0902.1089 [astroph.IM]].
 [7] M. Actis et al. [CTA Consortium Collaboration], Design concepts for the Cherenkov Telescope Array CTA: An advanced facility for groundbased highenergy gammaray astronomy, Exper. Astron. 32 (2011) 193 [arXiv:1008.3703 [astroph.IM]].
 [8] N. Bernal and S. PalomaresRuiz, Constraining Dark Matter Properties with GammaRays from the Galactic Center with FermiLAT, Nucl. Phys. B 857 (2012) 380 [arXiv:1006.0477 [astroph.HE]].
 [9] N. Bernal and S. PalomaresRuiz, Constraining the Milky Way Dark Matter Density Profile with GammaRays with FermiLAT, JCAP 1201 (2012) 006 [arXiv:1103.2377 [astroph.HE]].
 [10] K. Griest and D. Seckel, Three exceptions in the calculation of relic abundances, Phys. Rev. D 43 (1991) 3191. doi:10.1103/PhysRevD.43.3191
 [11] P. Gondolo and G. Gelmini, Cosmic abundances of stable particles: Improved analysis, Nucl. Phys. B 360 (1991) 145. doi:10.1016/05503213(91)904384
 [12] L. Roszkowski, E. M. Sessolo, S. Trojanowski and A. J. Williams, Reconstructing WIMP properties through an interplay of signal measurements in direct detection, FermiLAT, and CTA searches for dark matter, JCAP 1608 (2016) no.08, 033 [arXiv:1603.06519 [astroph.CO]].
 [13] M. L. Ahnen et al. [MAGIC and FermiLAT Collaborations], Limits to dark matter annihilation crosssection from a combined analysis of MAGIC and FermiLAT observations of dwarf satellite galaxies, JCAP 1602 (2016) no.02, 039 [arXiv:1601.06590 [astroph.HE]].
 [14] A. Albert et al. [FermiLAT and DES Collaborations], Searching for Dark Matter Annihilation in Recently Discovered Milky Way Satellites with FermiLAT, Astrophys. J. 834 (2017) no.2, 110 [arXiv:1611.03184 [astroph.HE]].
 [15] S. Profumo, F. S. Queiroz and C. E. Yaguna, Extending FermiLAT and H.E.S.S. Limits on Gammaray Lines from Dark Matter Annihilation, Mon. Not. Roy. Astron. Soc. 461 (2016) no.4, 3976 [arXiv:1602.08501 [astroph.HE]].
 [16] H. Abdallah et al. [HESS Collaboration], Search for dark matter annihilations towards the inner Galactic halo from 10 years of observations with H.E.S.S, Phys. Rev. Lett. 117 (2016) no.11, 111301 doi:10.1103/PhysRevLett.117.111301 [arXiv:1607.08142 [astroph.HE]].
 [17] J. Hisano, S. Matsumoto and M. M. Nojiri, Phys. Rev. D 67 (2003) 075014 doi:10.1103/PhysRevD.67.075014 [hepph/0212022].
 [18] J. Hisano, S. Matsumoto and M. M. Nojiri, Phys. Rev. Lett. 92 (2004) 031303 doi:10.1103/PhysRevLett.92.031303 [hepph/0307216].
 [19] J. Hisano, S. Matsumoto, M. M. Nojiri and O. Saito, Nonperturbative effect on dark matter annihilation and gamma ray signature from galactic center, Phys. Rev. D 71 (2005) 063528 [hepph/0412403].
 [20] M. Cirelli, A. Strumia and M. Tamburini, Nucl. Phys. B 787 (2007) 152 doi:10.1016/j.nuclphysb.2007.07.023 [arXiv:0706.4071 [hepph]].
 [21] N. ArkaniHamed, D. P. Finkbeiner, T. R. Slatyer and N. Weiner, Phys. Rev. D 79 (2009) 015014 doi:10.1103/PhysRevD.79.015014 [arXiv:0810.0713 [hepph]].
 [22] S. Cassel, J. Phys. G 37 (2010) 105009 [arXiv:0903.5307 [hepph]].
 [23] T. R. Slatyer, JCAP 1002 (2010) 028 [arXiv:0910.5713 [hepph]].
 [24] J. L. Feng, M. Kaplinghat and H. B. Yu, Phys. Rev. D 82 (2010) 083525 [arXiv:1005.4678 [hepph]].
 [25] S. El Hedri, A. Kaminska and M. de Vries, arXiv:1612.02825 [hepph].
 [26] K. Blum, R. Sato and T. R. Slatyer, JCAP 1606 (2016) no.06, 021 [arXiv:1603.01383 [hepph]].
 [27] J. Fan and M. Reece, In Wino Veritas? Indirect Searches Shed Light on Neutralino Dark Matter, JHEP 1310 (2013) 124 [arXiv:1307.4400 [hepph]].
 [28] T. Cohen, M. Lisanti, A. Pierce and T. R. Slatyer, Wino Dark Matter Under Siege, JCAP 1310 (2013) 061 [arXiv:1307.4082 [hepph]].
 [29] A. Hryczuk, I. Cholis, R. Iengo, M. Tavakoli and P. Ullio, Indirect Detection Analysis: Wino Dark Matter Case Study, JCAP 1407 (2014) 031 [arXiv:1401.6212 [astroph.HE]].
 [30] G. Duda, G. Gelmini, P. Gondolo, J. Edsjo and J. Silk, Indirect detection of a subdominant density component of cold dark matter, Phys. Rev. D 67 (2003) 023505 doi:10.1103/PhysRevD.67.023505 [hepph/0209266].
 [31] M. Beltran, D. Hooper, E. W. Kolb and Z. C. Krusberg, Deducing the nature of dark matter from direct and indirect detection experiments in the absence of collider signatures of new physics, Phys. Rev. D 80 (2009) 043509 [arXiv:0808.3384 [hepph]].
 [32] Q. H. Cao, C. R. Chen, C. S. Li and H. Zhang, Effective Dark Matter Model: Relic density, CDMS II, Fermi LAT and LHC, JHEP 1108 (2011) 018 [arXiv:0912.4511 [hepph]].
 [33] M. Beltran, D. Hooper, E. W. Kolb, Z. A. C. Krusberg and T. M. P. Tait, Maverick dark matter at colliders, JHEP 1009 (2010) 037 [arXiv:1002.4137 [hepph]].
 [34] J. Goodman, M. Ibe, A. Rajaraman, W. Shepherd, T. M. P. Tait and H. B. Yu, Constraints on Dark Matter from Colliders, Phys. Rev. D 82 (2010) 116010 [arXiv:1008.1783 [hepph]].
 [35] J. Goodman, M. Ibe, A. Rajaraman, W. Shepherd, T. M. P. Tait and H. B. Yu, Gamma Ray Line Constraints on Effective Theories of Dark Matter, Nucl. Phys. B 844 (2011) 55 [arXiv:1009.0008 [hepph]].
 [36] J. Kumar and D. Marfatia, Matrix element analyses of dark matter scattering and annihilation, Phys. Rev. D 88 (2013) no.1, 014035 [arXiv:1305.1611 [hepph]].
 [37] J. Goodman and W. Shepherd, LHC Bounds on UVComplete Models of Dark Matter, arXiv:1111.2359 [hepph].
 [38] M. T. Frandsen, F. Kahlhoefer, A. Preston, S. Sarkar and K. SchmidtHoberg, LHC and Tevatron Bounds on the Dark Matter Direct Detection CrossSection for Vector Mediators, JHEP 1207 (2012) 123 [arXiv:1204.3839 [hepph]].
 [39] H. Dreiner, D. Schmeier and J. Tattersall, Contact Interactions Probe Effective Dark Matter Models at the LHC, Europhys. Lett. 102 (2013) no.5, 51001 [arXiv:1303.3348 [hepph]].
 [40] O. Buchmueller, M. J. Dolan and C. McCabe, Beyond Effective Field Theory for Dark Matter Searches at the LHC, JHEP 1401 (2014) 025 [arXiv:1308.6799 [hepph]].
 [41] M. Abdullah, A. DiFranzo, A. Rajaraman, T. M. P. Tait, P. Tanedo and A. M. Wijangco, Hidden onshell mediators for the Galactic Center ray excess, Phys. Rev. D 90 (2014) 035004 [arXiv:1404.6528 [hepph]].
 [42] C. Karwin, S. Murgia, T. M. P. Tait, T. A. Porter and P. Tanedo, Dark Matter Interpretation of the FermiLAT Observation Toward the Galactic Center, arXiv:1612.05687 [hepph].
 [43] F. D’Eramo, B. J. Kavanagh and P. Panci, Probing Leptophilic Dark Sectors with Hadronic Processes, arXiv:1702.00016 [hepph].
 [44] G. Aad et al. [ATLAS Collaboration], Search for dark matter candidates and large extra dimensions in events with a jet and missing transverse momentum with the ATLAS detector, JHEP 1304 (2013) 075 [arXiv:1210.4491 [hepex]].
 [45] V. Khachatryan et al. [CMS Collaboration], Search for dark matter, extra dimensions, and unparticles in monojet events in protonâproton collisions at TeV, Eur. Phys. J. C 75 (2015) no.5, 235 [arXiv:1408.3583 [hepex]].
 [46] J. Fan, M. Reece and L. T. Wang, Nonrelativistic effective theory of dark matter direct detection, JCAP 1011 (2010) 042 [arXiv:1008.1591 [hepph]].
 [47] A. L. Fitzpatrick, W. Haxton, E. Katz, N. Lubbers and Y. Xu, The Effective Field Theory of Dark Matter Direct Detection, JCAP 1302 (2013) 004 [arXiv:1203.3542 [hepph]].
 [48] F. D’Eramo and M. Procura, Connecting Dark Matter UV Complete Models to Direct Detection Rates via Effective Field Theory, JHEP 1504 (2015) 054 [arXiv:1411.3342 [hepph]].
 [49] M. Srednicki, R. Watkins and K. A. Olive, Calculations of Relic Densities in the Early Universe, Nucl. Phys. B 310 (1988) 693.
 [50] J. M. Zheng, Z. H. Yu, J. W. Shao, X. J. Bi, Z. Li and H. H. Zhang, Constraining the interaction strength between dark matter and visible matter: I. fermionic dark matter, Nucl. Phys. B 854 (2012) 350 [arXiv:1012.2022 [hepph]].
 [51] F. D’Eramo, B. J. Kavanagh and P. Panci, You can hide but you have to run: direct detection with vector mediators, JHEP 1608 (2016) 111 [arXiv:1605.04917 [hepph]].
 [52] D. S. Akerib et al. [LUX Collaboration], Results from a search for dark matter in the complete LUX exposure, Phys. Rev. Lett. 118 (2017) no.2, 021303 [arXiv:1608.07648 [astroph.CO]].
 [53] H. P. Nilles and S. Raby, Supersymmetry and the strong CP problem, Nucl. Phys. B 198 (1982) 102.
 [54] J. M. Frere and J. M. Gerard, Axions and Supersymmetry, Lett. Nuovo Cim. 37 (1983) 135.
 [55] K. Tamvakis and D. Wyler, Broken Global Symmetries in Supersymmetric Theories, Phys. Lett. B 112 (1982) 451.
 [56] L. Covi, J. E. Kim and L. Roszkowski, Axinos as cold dark matter, Phys. Rev. Lett. 82 (1999) 4180 [hepph/9905212].
 [57] L. Covi, H. B. Kim, J. E. Kim and L. Roszkowski, Axinos as dark matter, JHEP 0105 (2001) 033 [hepph/0101009].
 [58] J. E. Kim and M. S. Seo, Mixing of axino and goldstino, and axino mass, Nucl. Phys. B 864 (2012) 296 [arXiv:1204.5495 [hepph]].
 [59] K. Y. Choi, J. E. Kim and L. Roszkowski, Review of axino dark matter, J. Korean Phys. Soc. 63 (2013) 1685 [arXiv:1307.3330 [astroph.CO]].
 [60] S. Deser and B. Zumino, Broken Supersymmetry and Supergravity, Phys. Rev. Lett. 38 (1977) 1433.
 [61] E. Cremmer, B. Julia, J. Scherk, P. van Nieuwenhuizen, S. Ferrara and L. Girardello, Superhiggs effect in supergravity with general scalar interactions, Phys. Lett. B 79 (1978) 231.
 [62] E. Cremmer, B. Julia, J. Scherk, S. Ferrara, L. Girardello and P. van Nieuwenhuizen, Spontaneous Symmetry Breaking and Higgs Effect in Supergravity Without Cosmological Constant, Nucl. Phys. B 147 (1979) 105.
 [63] E. Cremmer, S. Ferrara, L. Girardello and A. Van Proeyen, YangMills Theories with Local Supersymmetry: Lagrangian, Transformation Laws and SuperHiggs Effect, Nucl. Phys. B 212 (1983) 413.
 [64] J. L. Feng, A. Rajaraman and F. Takayama, Superweakly interacting massive particles, Phys. Rev. Lett. 91 (2003) 011302 [hepph/0302215].
 [65] M. Bolz, A. Brandenburg and W. Buchmuller, Thermal production of gravitinos, Nucl. Phys. B 606 (2001) 518 Erratum: [Nucl. Phys. B 790 (2008) 336] [hepph/0012052].
 [66] A. Brandenburg and F. D. Steffen, Axino dark matter from thermal production, JCAP 0408 (2004) 008 [hepph/0405158].
 [67] J. Pradler and F. D. Steffen, Thermal gravitino production and collider tests of leptogenesis, Phys. Rev. D 75 (2007) 023509 [hepph/0608344].
 [68] V. S. Rychkov and A. Strumia, Thermal production of gravitinos, Phys. Rev. D 75 (2007) 075011 [hepph/0701104].
 [69] K. J. Bae, K. Choi and S. H. Im, Effective Interactions of Axion Supermultiplet and Thermal Production of Axino Dark Matter, JHEP 1108 (2011) 065 [arXiv:1106.2452 [hepph]].
 [70] K. Y. Choi, L. Covi, J. E. Kim and L. Roszkowski, Axino Cold Dark Matter Revisited, JHEP 1204 (2012) 106 [arXiv:1108.2282 [hepph]].
 [71] K. J. Bae, E. J. Chun and S. H. Im, Cosmology of the DFSZ axino, JCAP 1203 (2012) 013 [arXiv:1111.5962 [hepph]].
 [72] M. H. Reno and D. Seckel, Primordial Nucleosynthesis: The Effects of Injecting Hadrons, Phys. Rev. D 37 (1988) 3441.
 [73] M. Kawasaki and T. Moroi, Electromagnetic cascade in the early universe and its application to the big bang nucleosynthesis, Astrophys. J. 452 (1995) 506 [astroph/9412055].
 [74] J. L. Feng, A. Rajaraman and F. Takayama, SuperWIMP dark matter signals from the early universe, Phys. Rev. D 68 (2003) 063504 [hepph/0306024].
 [75] J. L. Feng, S. f. Su and F. Takayama, SuperWIMP gravitino dark matter from slepton and sneutrino decays, Phys. Rev. D 70 (2004) 063514 [hepph/0404198].
 [76] K. Jedamzik, M. Lemoine and G. Moultaka, Gravitino, axino, KaluzaKlein graviton warm and mixed dark matter and reionisation, JCAP 0607 (2006) 010 doi:10.1088/14757516/2006/07/010 [astroph/0508141].
 [77] J. R. Ellis, J. E. Kim and D. V. Nanopoulos, Cosmological Gravitino Regeneration and Decay, Phys. Lett. B 145 (1984) 181.
 [78] W. Hu and J. Silk, Thermalization constraints and spectral distortions for massive unstable relic particles, Phys. Rev. Lett. 70 (1993) 2661.
 [79] S. Ambrosanio, G. L. Kane, G. D. Kribs, S. P. Martin and S. Mrenna, Search for supersymmetry with a light gravitino at the Fermilab Tevatron and CERN LEP colliders, Phys. Rev. D 54 (1996) 5395 [hepph/9605398].
 [80] W. Buchmuller, K. Hamaguchi, M. Ratz and T. Yanagida, Supergravity at colliders, Phys. Lett. B 588 (2004) 90 [hepph/0402179].
 [81] L. Covi, L. Roszkowski, R. Ruiz de Austri and M. Small, Axino dark matter and the CMSSM, JHEP 0406 (2004) 003 [hepph/0402240].
 [82] J. L. Feng, S. Su and F. Takayama, Supergravity with a gravitino LSP, Phys. Rev. D 70 (2004) 075019 [hepph/0404231].
 [83] L. Roszkowski, R. Ruiz de Austri and K. Y. Choi, Gravitino dark matter in the CMSSM and implications for leptogenesis and the LHC, JHEP 0508 (2005) 080 [hepph/0408227].
 [84] K. Hamaguchi, Y. Kuno, T. Nakaya and M. M. Nojiri, A Study of late decaying charged particles at future colliders, Phys. Rev. D 70 (2004) 115007 [hepph/0409248].
 [85] J. L. Feng and B. T. Smith, Slepton trapping at the large hadron and international linear colliders, Phys. Rev. D 71 (2005) 015004 Erratum: [Phys. Rev. D 71 (2005) 019904] [hepph/0409278].
 [86] A. Brandenburg, L. Covi, K. Hamaguchi, L. Roszkowski and F. D. Steffen, Signatures of axinos and gravitinos at colliders, Phys. Lett. B 617 (2005) 99 [hepph/0501287].
 [87] K. Y. Choi, L. Roszkowski and R. Ruiz de Austri, Determining reheating temperature at colliders with axino or gravitino dark matter, JHEP 0804 (2008) 016 [arXiv:0710.3349 [hepph]].
 [88] J. L. Feng, M. Kamionkowski and S. K. Lee, Light Gravitinos at Colliders and Implications for Cosmology, Phys. Rev. D 82 (2010) 015012 [arXiv:1004.4213 [hepph]].
 [89] A. Arbey, M. Battaglia, L. Covi, J. Hasenkamp and F. Mahmoudi, LHC constraints on Gravitino Dark Matter, Phys. Rev. D 92 (2015) no.11, 115008 [arXiv:1505.04595 [hepph]].
 [90] M. Ibe, M. Suzuki and T. T. Yanagida, Revisiting gravitino dark matter in thermal leptogenesis, JHEP 1702 (2017) 063 [arXiv:1609.06834 [hepph]].
 [91] G. Aad et al. [ATLAS Collaboration], Search for nonpointing and delayed photons in the diphoton and missing transverse momentum final state in 8 TeV collisions at the LHC using the ATLAS detector, Phys. Rev. D 90 (2014) no.11, 112005 [arXiv:1409.5542 [hepex]].
 [92] F. Maltoni, A. Martini, K. Mawatari and B. Oexl, Signals of a superlight gravitino at the LHC, JHEP 1504 (2015) 021 [arXiv:1502.01637 [hepph]].
 [93] G. Aad et al. [ATLAS Collaboration], Search for new phenomena in final states with an energetic jet and large missing transverse momentum in pp collisions at 8 TeV with the ATLAS detector, Eur. Phys. J. C 75 (2015) no.7, 299 Erratum: [Eur. Phys. J. C 75 (2015) no.9, 408] [arXiv:1502.01518 [hepex]].
 [94] V. Khachatryan et al. [CMS Collaboration], Search for exotic decays of a Higgs boson into undetectable particles and one or more photons, Phys. Lett. B 753 (2016) 363 [arXiv:1507.00359 [hepex]].
 [95] J. C. Park, S. C. Park and K. Kong, Xray line signal from 7 keV axino dark matter decay, Phys. Lett. B 733 (2014) 217 [arXiv:1403.1536 [hepph]].
 [96] K. Y. Choi and O. Seto, Xray line signal from decaying axino warm dark matter, Phys. Lett. B 735 (2014) 92 [arXiv:1403.1782 [hepph]].
 [97] S. P. Liew, Axino dark matter in light of an anomalous Xray line, JCAP 1405 (2014) 044 [arXiv:1403.6621 [hepph], arXiv:1403.6621].
 [98] N.E. Bomark and L. Roszkowski, 3.5 keV xray line from decaying gravitino dark matter, Phys. Rev. D 90 (2014) 011701 [arXiv:1403.6503 [hepph]].
 [99] G. Arcadi, L. Covi and M. Nardecchia, Gravitino Dark Matter and lowscale Baryogenesis, Phys. Rev. D 92 (2015) no.11, 115006 [arXiv:1507.05584 [hepph]].
 [100] G. A. GomezVargas, D. E. LopezFogliani, C. Munoz, A. D. Perez and R. Ruiz de Austri, Search for sharp and smooth spectral signatures of SSM gravitino dark matter with FermiLAT, arXiv:1608.08640 [hepph].
 [101] L. Roszkowski, S. Trojanowski and K. Turzynski, Axino dark matter with low reheating temperature, JHEP 1511 (2015) 139 [arXiv:1507.06164 [hepph]].
 [102] G. F. Giudice, E. W. Kolb and A. Riotto, Largest temperature of the radiation era and its cosmological implications, Phys. Rev. D 64 (2001) 023508 [hepph/0005123].
 [103] L. J. Hall, K. Jedamzik, J. MarchRussell and S. M. West, FreezeIn Production of FIMP Dark Matter, JHEP 1003 (2010) 080 [arXiv:0911.1120 [hepph]].
 [104] L. Roszkowski, S. Trojanowski and K. Turzyński, Neutralino and gravitino dark matter with low reheating temperature, JHEP 1411 (2014) 146 [arXiv:1406.0012 [hepph]].
 [105] J. E. Kim, Weak Interaction Singlet and Strong CP Invariance, Phys. Rev. Lett. 43 (1979) 103.
 [106] M. A. Shifman, A. I. Vainshtein and V. I. Zakharov, Can Confinement Ensure Natural CP Invariance of Strong Interactions?, Nucl. Phys. B 166 (1980) 493.
 [107] K. Y. Choi, J. E. Kim, H. M. Lee and O. Seto, Neutralino dark matter from heavy axino decay, Phys. Rev. D 77 (2008) 123501 [arXiv:0801.0491 [hepph]].
 [108] K. J. Bae, J. H. Huh and J. E. Kim, Update of axion CDM energy, JCAP 0809 (2008) 005 [arXiv:0806.0497 [hepph]].
 [109] M. Kawasaki, K. Kohri and T. Moroi, BigBang nucleosynthesis and hadronic decay of longlived massive particles, Phys. Rev. D 71 (2005) 083502 [astroph/0408426].
 [110] K. Jedamzik, Big bang nucleosynthesis constraints on hadronically and electromagnetically decaying relic neutral particles, Phys. Rev. D 74 (2006) 103509 [hepph/0604251].
 [111] M. Ibe, A. Kamada and S. Matsumoto, Imprints of nonthermal Wino dark matter on smallscale structure, Phys. Rev. D 87 (2013) no.6, 063511 [arXiv:1210.0191 [hepph]].
 [112] M. Kawasaki, K. Kohri, T. Moroi and A. Yotsuyanagi, BigBang Nucleosynthesis and Gravitino, Phys. Rev. D 78 (2008) 065011 [arXiv:0804.3745 [hepph]].
 [113] J. F. Navarro, C. S. Frenk and S. D. M. White, The Structure of cold dark matter halos, Astrophys. J. 462 (1996) 563 [astroph/9508025].
 [114] D. Merritt, J. F. Navarro, A. Ludlow and A. Jenkins, A Universal density profile for dark and luminous matter?, Astrophys. J. 624 (2005) L85 [astroph/0502515].
 [115] A. W. Graham, D. Merritt, B. Moore, J. Diemand and B. Terzic, Empirical Models for Dark Matter Halos. II. Inner profile slopes, dynamical profiles, and , Astron. J. 132 (2006) 2701 [astroph/0608613].
 [116] A. Burkert, The Structure of dark matter halos in dwarf galaxies, IAU Symp. 171 (1996) 175 [Astrophys. J. 447 (1995) L25] [astroph/9504041].
 [117] M. Cirelli et al., PPPC 4 DM ID: A Poor Particle Physicist Cookbook for Dark Matter Indirect Detection, JCAP 1103 (2011) 051 Erratum: [JCAP 1210 (2012) E01] [arXiv:1012.4515 [hepph]].
 [118] V. Lefranc, G. A. Mamon and P. Panci, Prospects for annihilating Dark Matter towards Milky Way’s dwarf galaxies by the Cherenkov Telescope Array, JCAP 1609 (2016) no.09, 021 [arXiv:1605.02793 [astroph.HE]].
 [119] H. Zhao, Analytical models for galactic nuclei, Mon. Not. Roy. Astron. Soc. 278 (1996) 488 [astroph/9509122].
 [120] M. Ackermann et al. [FermiLAT Collaboration], Dark matter constraints from observations of 25 Milky Way satellite galaxies with the Fermi Large Area Telescope, Phys. Rev. D 89 (2014) 042001 [arXiv:1310.0828 [astroph.HE]].
 [121] T. MarrodÃ¡n Undagoitia and L. Rauch, Dark matter directdetection experiments, J. Phys. G 43 (2016) no.1, 013001 [arXiv:1509.08767 [physics.insdet]].
 [122] J. R. Ellis, R. A. Flores and J. D. Lewin, Rates for Inelastic Nuclear Excitation by Dark Matter Particles, Phys. Lett. B 212 (1988) 375.
 [123] J. Engel and P. Vogel, Neutralino inelastic scattering with subsequent detection of nuclear gammarays, Phys. Rev. D 61 (2000) 063503 [hepph/9910409].
 [124] M. Pato, What can(not) be measured with tonscale dark matter direct detection experiments, JCAP 1110 (2011) 035 [arXiv:1106.0743 [astroph.CO]].
 [125] J. L. Newstead, T. D. Jacques, L. M. Krauss, J. B. Dent and F. Ferrer, Scientific reach of multitonscale dark matter direct detection experiments, Phys. Rev. D 88 (2013) no.7, 076011 [arXiv:1306.3244 [astroph.CO]].
 [126] A. H. G. Peter, V. Gluscevic, A. M. Green, B. J. Kavanagh and S. K. Lee, WIMP physics with ensembles of directdetection experiments, Phys. Dark Univ. 56 (2014) 45 [arXiv:1310.7039 [astroph.CO]].
 [127] N. Bernal, A. Goudelis, Y. Mambrini and C. Munoz, Determining the WIMP mass using the complementarity between direct and indirect searches and the ILC, JCAP 0901 (2009) 046 [arXiv:0804.1976 [hepph]].
 [128] C. Arina, G. Bertone and H. Silverwood, Complementarity of direct and indirect Dark Matter detection experiments, Phys. Rev. D 88 (2013) no.1, 013002 [arXiv:1304.5119 [hepph]].
 [129] B. J. Kavanagh, M. Fornasa and A. M. Green, Probing WIMP particle physics and astrophysics with direct detection and neutrino telescope data, Phys. Rev. D 91 (2015) no.10, 103533 [arXiv:1410.8051 [astroph.CO]].
 [130] P. Ciafaloni, D. Comelli, A. Riotto, F. Sala, A. Strumia and A. Urbano, Weak Corrections are Relevant for Dark Matter Indirect Detection, JCAP 1103 (2011) 019 [arXiv:1009.0224 [hepph]].
 [131] J. Buch, M. Cirelli, G. Giesen and M. Taoso, PPPC 4 DM secondary: A Poor Particle Physicist Cookbook for secondary radiation from Dark Matter, JCAP 1509 (2015) no.09, 037 [arXiv:1505.01049 [hepph]].
 [132] G. Blanger, F. Boudjema, A. Pukhov and A. Semenov, micrOMEGAs4.1: two dark matter candidates, Comput. Phys. Commun. 192 (2015) 322 [arXiv:1407.6129 [hepph]].
 [133] A. Djouadi, M. M. Muhlleitner and M. Spira, Decays of supersymmetric particles: The Program SUSYHIT (SUspectSdecaYHdecayInTerface), Acta Phys. Polon. B 38 (2007) 635 [hepph/0609292].
 [134] J. Hisano, S. Matsumoto, M. Nagai, O. Saito and M. Senami, Nonperturbative effect on thermal relic abundance of dark matter, Phys. Lett. B 646 (2007) 34 [hepph/0610249].
 [135] A. Hryczuk and R. Iengo, The oneloop and Sommerfeld electroweak corrections to the Wino dark matter annihilation, JHEP 1201 (2012) 163 Erratum: [JHEP 1206 (2012) 137] [arXiv:1111.2916 [hepph]].
 [136] M. Beneke, A. Bharucha, F. Dighera, C. Hellmann, A. Hryczuk, S. Recksiegel and P. RuizFemenia, Relic density of winolike dark matter in the MSSM, JHEP 1603 (2016) 119 [arXiv:1601.04718 [hepph]].
 [137] G. Ovanesyan, N. L. Rodd, T. R. Slatyer and I. W. Stewart, Oneloop correction to heavy dark matter annihilation, Phys. Rev. D 95 (2017) no.5, 055001 [arXiv:1612.04814 [hepph]].
 [138] A. Hryczuk, The Sommerfeld enhancement for scalar particles and application to sfermion coannihilation regions, Phys. Lett. B 699 (2011) 271 [arXiv:1102.4295 [hepph]].
 [139] G. Ovanesyan, T. R. Slatyer and I. W. Stewart, Heavy Dark Matter Annihilation from Effective Field Theory, Phys. Rev. Lett. 114 (2015) no.21, 211302 [arXiv:1409.8294 [hepph]].
 [140] J. Hisano, S. Matsumoto, M. M. Nojiri and O. Saito, Direct detection of the Wino and Higgsinolike neutralino dark matters at oneloop level, Phys. Rev. D 71 (2005) 015007 [hepph/0407168].
 [141] F. Feroz, M. Hobson, M. Bridges, MultiNest: an efficient and robust Bayesian inference tool for cosmology and particle physics, Mon. Not. Roy. Astron. Soc. 398 (2009) 16011614 arXiv:0809.3437 [astroph].
 [142] B. Allanach, SOFTSUSY: a program for calculating supersymmetric spectra , Comput. Phys. Commun. 143 (2002) 305331 [hepph/0104145].
 [143] G. Aad et al. [ATLAS and CMS Collaborations], Combined Measurement of the Higgs Boson Mass in Collisions at and 8 TeV with the ATLAS and CMS Experiments, Phys. Rev. Lett. 114 (2015) 191803 [arXiv:1503.07589 [hepex]].
 [144] K. Kowalska, Phenomenological MSSM in light of new 13 TeV LHC data, Eur. Phys. J. C 76 (2016) no.12, 684 [arXiv:1608.02489 [hepph]].
 [145] A. Arbey, F. Mahmoudi, SuperIso Relic: A program for calculating relic density and flavor physics observables in Supersymmetry, Comput. Phys. Commun. 181 (2010) 12771292 arXiv:0906.0369 [hepph].
 [146] C. Patrignani, Review of Particle Physics, Chin. Phys. C 40 (2016) no.10, 100001.
 [147] Y. Amhis et al. [Heavy Flavor Averaging Group (HFAG) Collaboration], Averages of hadron, hadron, and lepton properties as of summer 2014, arXiv:1412.7515 [hepex].
 [148] http://www.slac.stanford.edu/xorg/hfag/rare/2012/radll/index.ht
 [149] T. Aaltonen et al. [CDF Collaboration], Search for and decays with the full CDF Run II data set, Phys. Rev. D 87 (2013) no.7, 072003 [arXiv:1301.7048 [hepex]].
 [150] V. Khachatryan et al. [CMS and LHCb Collaborations], Observation of the rare decay from the combined analysis of CMS and LHCb data, Nature 522 (2015) 68 [arXiv:1411.4413 [hepex]].
 [151] E. J. Chun, Dark matter in the KimNilles mechanism, Phys. Rev. D 84 (2011) 043509 [arXiv:1104.2219 [hepph]].