#
-wave Annihilating Dark Matter from a Decaying Predecessor

and the Galactic Center Excess

###### Abstract

Dark matter (DM) annihilations have been widely studied as a possible explanation of excess gamma rays from the galactic center seen by Fermi/LAT. However most such models are in conflict with constraints from dwarf spheroidals. Motivated by this tension, we show that -wave annihilating dark matter can easily accommodate both sets of observations due to the lower DM velocity dispersion in dwarf galaxies. Explaining the DM relic abundance is then challenging. We outline a scenario in which the usual thermal abundance is obtained through -wave annihilations of a metastable particle, that eventually decays into the -wave annihilating DM of the present epoch. The couplings and lifetime of the decaying particle are constrained by big bang nucleosynthesis, the cosmic microwave background and direct detection, but significant regions of parameter space are viable. A sufficiently large -wave cross section can be found by annihilation into light mediators, that also give rise to Sommerfeld enhancement. A prediction of the scenario is enhanced annihilations in galaxy clusters.

Fermi-LAT Observations of the galactic center (GC) provide evidence of a gamma-ray excess in the multi-GeV energy range Goodenough:2009gk (); Hooper:2010mq (); Hooper:2011ti (); Abazajian:2012pn (); Zhou:2014lva (); Calore:2014xka (); Daylan:2014rsa (); FermiData (); TheFermi-LAT:2015kwa (). Millisecond pulsars (MSPs) are a favored astrophysical source to explain the signal Lee:2015fea (); Bartels:2015aea (); O'Leary:2015gfa (), but there is debate in the literature as to whether the required numbers of MSPs in the GC for this purpose is consistent with the number resolved by Fermi Linden:2015qha () or expected on theoretical or empirical grounds Cholis:2014lta (). Other possible astrophysical explanations have been presented Carlson:2014cwa (); Petrovic:2014uda (); Cholis:2015dea (); Gaggero:2015nsa (); deBoer:2015kta (); Carlson:2016iis (), but dark matter annihilation into charged particles that lead to gamma rays remains a possibility that has attracted great interest. Further data should eventually be able to distinguish between the different possibilities Yuan:2014yda (); Lee:2014mza ().

There is tension between most dark matter (DM) explanations of the
galactic center excess (GCE) and constraints on dark matter annihilation coming from
observations of dwarf spheroidal galaxies Ackermann:2015zua (); Abazajian:2015raa ().^{4}^{4}4Analyses of known dwarf spheroidal galaxies have revealed no significant excess gamma-ray emission. However, there have been claims of possible signals from the recently discovered Bechtol:2015cbp (); Drlica-Wagner:2015ufc () dwarf spheroidal candidates Reticulum II Geringer-Sameth:2015lua (); Hooper:2015ula () and Turcana III Li:2015kag (). These results are somewhat in dispute, with a Fermi-LAT analysis of Reticulum II using more data Drlica-Wagner:2015xua () claiming no significant excess.
(Further complementary constraints come from searches for GeV emission
in the large Magellanic cloud Buckley:2015doa () or subhalos
of the Milky Way
Bertoni:2015mla ().)
The best fits for DM mass and annihilation cross section for the
GCE lie in regions that tend to be excluded by factors of a
few by the dwarf spheroidal limits. A possible way of alleviating
this tension is to assume that the annihilation is into electrons,
a scenario in which the GCE is primarily produced through inverse Compton scattering which is suppressed in dwarfs because of their dilute
radiation fields
Liu:2014cma (); Calore:2014nla (); Kaplinghat:2015gha (). An additional idea to explain the discrepancy is a model of asymmetric DM where anti-DM is produced at late times via decays, leading to particles with enough kinetic energy to escape a dwarf galaxy but not the galactic center, where they annihilate with DM particles Hardy:2014dea ().

In this work we explore a different possibility, noting that the tension can be avoided if the dark matter annihilation rates are velocity-dependent. Since the velocity dispersion in the galactic center is significantly higher than that in dwarf galaxies, the GCE can be consistent with the lack of signals from dwarf spheroidal galaxies provided that the annihilation cross section increases with velocity. This is the case in models where -wave annihilations dominate, which is the subject of the present work. This scenario has recently been explored Zhao:2016xie () to alleviate tension between the dwarf spheroidal constraints and DM explanations of the AMS-02 positron excess. We take a similar approach for the GCE. An immediate challenge is how to obtain the right relic density since the cross section needed for the GCE is of order , the usual value associated with a thermal origin for the relic density. But if has such a value in the GC today, it would have been orders of magnitude larger in the early universe, leading to a negligible thermal abundance. We address this by showing how the current generation of -wave annihilating dark matter could have arisen through the decays of a metastable predecessor DM particle that has a thermal origin. The decays can take place at temperatures ranging from eV to several GeV. By this time the -wave annihilations would be out of equilibrium despite their relatively large cross section.

The annihilation cross section needed to explain the GCE requires large couplings to compensate for the -wave suppresssion. Such large copulings would generically tend to also give strong interactions of dark matter with nuclei. However constraints from direct detection can be satisfied if the dark matter annihilates into light mediators Abdullah:2014lla (); Elor:2015tva () that subsequently decay into standard model particles. The light mediators also lead to Sommerfeld enhanced annihilation, allowing us to avoid nonperturbatively large couplings. In this way we are able to find viable models that have reasonably small couplings.

In section I we parametrize the -wave annihilation cross section in the Milky Way (MW) and in dwarf spheroidal galaxies, in terms of assumed velocity dispersion profiles, leading to modified -factors that are relevant for comparison to observations. In section II we give the results of the galactic propagation simulations used to compute the expected signal from the galactic center, including the effects of inverse Compton scattering and Bremsstrahlung radiation. This yields fits to the data in the plane of DM mass versus annihilation cross section . We then derive upper limits on in the same plane from dwarf spheroidals and galaxy clusters. In section III we show that -wave annihilations of the desired strength would lead to strong suppression of the DM abundance at freeze-out, unless some nonthermal origin prevails. Here we present the scenario of decaying DM whose density is determined by the usual -wave process, and the conditions under which this provides a consistent description. Three examples of decay channels leading to different phenomenology are presented, to illustrate the range of possibilities. In section IV we systematically explore observational constraints on these models coming from cosmology, astrophysical line searches, direct searches, and colliders. In section V we provide a concrete model of annihilation into light scalar mediators to show that the desired large cross section can be achieved with reasonable values of the couplings in a renormalizable model. Conclusions are given in section VI.

## I Annihilation Cross Section

The expected signal from either the GC or dwarf spheroidals is proportional to the phase-space averaged cross section,

(1) |

for a velocity distribution , where is the relative velocity between the two annihilating particles and the escape velocity depends upon radial position in the galaxy. In this work we consider Dirac fermion dark matter. Self-annihilating Majorana dark matter would introduce an additional factor of into equation (1). Following Robertson:2009bh () and others, we adopt a Maxwell-Boltzmann distribution,

(2) |

where is the velocity dispersion at the given . The normalization factor in (2) is appropriate in the limit of large escape velocity, . Numerically we find that this approximation is well-suited to the present applications.

We will be interested in -wave annihilation for which at low velocities with a constant. The phase-space averaged value is then

(3) |

In general is a function of . This dependence is potentially significant in the MW, unlike in dwarf spheroidals, whose radial dependence has been observed to be roughly constant. Regardless of these details, it is however clear that is several orders of magnitude lower in dwarf spheroidals (dSph) than in the MW if the cross section is -wave suppressed. Measured values of are less than 15 km/s in MW dSph satellites Walker:2007ju (), whereas most estimates of near the GC are km/s (see for example refs. Battaglia:2005rj (); Dehnen:2006cm ()). On the other hand, Fermi upper limits on from dSph observations are at most a factor of a few more stringent than the values of needed to fit the GCE.

### i.1 The Milky Way

The Milky Way, though composed predominantly of dark matter, has inner regions such as the bulge and bar (as well as Sagittarius A) which are dominated by baryonic matter or otherwise do not follow an Navarro-Frenk-White (NFW) Navarro:1995iw () profile. The velocity dispersion of dark matter in the Milky Way is difficult to measure directly in the inner region, hence we rely on simulations and theoretical estimates. In order to quantify the uncertainties associated with choosing a velocity dispersion profile, we base our profiles on the results of simulations RomanoDiaz:2008wz () that include baryonic matter to study the evolution of the Milky Way’s profile.

If the Milky Way contained no baryonic matter, it could be suitably modeled by an NFW profile. The resulting velocity dispersion, from fits to the aforementioned simulation, is

(4) |

with Chen:2009av (). When baryons are included, however, a slope of provides a better fit to the simulations Cirelli:2010nh (). We use the value , consistent with the results of Battaglia:2005rj (); Dehnen:2006cm ().

A second possibility that we consider is that the velocity dispersion of the Milky Way scales as a simple power law,

(5) |

as suggested by the results of ref. RomanoDiaz:2008wz (). A numerical fit to those results gives Chen:2009av () and, using our convention, a value of , which results in the same velocity dispersion at as eq. (4). Ref. RomanoDiaz:2008wz () resolves only down to radii kpc, so (5) need not hold at smaller radii. Nevertheless we extrapolate it to kpc to estimate the theoretical upper bound on the predicted GCE signal, which is greater for the ansatz (5) than for eq. (4). Since the observed signal is averaged over volume with weighting, the difference for the predicted GCE excess between the two assumptions is relatively small despite the fact that has very different behavior between the two as .

### i.2 Dwarf Spheroidal Galaxies

Dwarf spheroidal galaxies tend to have relatively flat observed velocity dispersion profiles out to large radii Walker:2007ju (). We therefore approximate them as being constant, independent of radius. In this case, the -factor for -wave annihilation is simply proportional to that for -wave. We define the former to be

(6) |

In ref. Ackermann:2013yva (). the -wave -factors of the 18 dwarf spheroidal galaxies for which kinematic data was available were computed. We use these to determine through the relation . Table 1 shows the velocity dispersions and -factors of the dwarf galaxies used.

Galaxy | (kms) | Ref. | ||
---|---|---|---|---|

Carina | 7.5 | 8.9 | Wilkinson:2006qq () | |

Draco | 13 | 10.1 | Wilkinson:2004fz () | |

Fornax | 11.1 | 9.3 | Walker:2005nt () | |

Leo I | 9.9 | 8.7 | Koch:2006in () | |

Leo II | 6.8 | 8.3 | Coleman:2007xe () | |

Sculptor | 9 | 9.6 | Westfall:2005ji () | |

Sextans | 8 | 9.3 | Wilkinson:2006qq () | |

Ursa Minor | 12 | 10.0 | Wilkinson:2004fz () | |

Bootes I | 6.6 | 9.5 | Belokurov:2006hf () | |

Canes Venatici I | 7.6 | 8.5 | Simon:2007dq () | |

Canes Venatcici II | 4.6 | 8.3 | Simon:2007dq () | |

Coma Berenices | 4.6 | 9.4 | Simon:2007dq () | |

Hercules | 5.1 | 8.6 | Simon:2007dq () | |

Leo IV | 3.3 | 8.0 | Simon:2007dq () | |

Segue 1 | 4.3 | 9.8 | Geha:2008zr () | |

Ursa Major I | 7.6 | 9.1 | Simon:2007dq () | |

Ursa Major II | 6.7 | 10.0 | Simon:2007dq () | |

Willman 1 | 4.0 | 9.3 | Willman:2010gy () |

## Ii Simulations and Indirect Limits

The observed gamma ray excess, if it originates from dark matter, can be the result of annihilations to SM particles. It has been shown that the observed flux can be fit by annihilations with a large branching ratio to , as would be expected for Higgs portal dark matter Daylan:2014rsa (); Calore:2014xka (); Cuoco:2016jqt (). Although most of these gamma rays are prompt (decay products of the quarks), a significant fraction comes from inverse Compton scattering and, to a lesser extent, from bremsstrahlung. While the prompt signal can be relatively easily computed, the ICS and bremsstrahlung contributions are more involved. To this end, we use the DRAGON Evoli:2008dv () code to simulate cosmic ray production and propagation from dark matter annihilations, and the GammaSky program which implements GALPROP Strong:1998fr () to simulate the ICS and bremsstrahlung contributions along the line of sight. GammaSky is as yet unreleased, but some results have been given in DiBernardo:2012zu ().

We have modified DRAGON to account for -wave annihilating dark matter, replacing the constant cross section appearing with the DM density by . We also incorporate a generalized NFW profile

(7) |

and the galactic diffusion parameters and magnetic field model used in ref. Calore:2014xka (), corresponding to their best-fit model (therein called Model F). The NFW parameters are taken to be (giving a local DM density of ), , and (the best-fit value for the GCE found in Calore:2014xka (); Daylan:2014rsa ()). The electron injection spectrum is taken from PPPC 4 Cirelli:2010xx (); Ciafaloni:2010ti (), as is the photon spectrum used in calculating the prompt contribution.

We will focus on models in which DM annihilates into on-shell scalar mediators that subsequently decay into SM particles, primarily . The prompt photon and electron spectra must be boosted with respect to those from DM annihilating at rest, to account for the velocity of when it decays. The decay spectrum into particles of type in the rest frame of the is denoted by . It is related to the spectrum in the center of mass frame of the system by Agrawal:2014oha (); Cline:2015qha ()

(8) |

where . This expression assumes that the final state particles are massless, which is approximately true for the electrons as well as the photons injected from decays.

The prompt photon spectrum can be calculated independently of the DRAGON simulation. Its integrated spectral flux (in units of ) is given by

(9) |

with defined in eq. (6). The total observed spectrum is equal to the sum of and the ICS+Bremsstrahlung spectrum determined from the simulations.

### ii.1 Simulation Results

We simulated the gamma ray flux for a range of dark matter masses () and compared the results to the GCE signals estimated in refs. Calore:2014xka (); Daylan:2014rsa (); FermiData (). The best fit regions are presented in fig. 1, which show the confidence intervals in the - plane for four different values for the mediator mass, GeV. The contours are generated by minimizing the of our simulated spectrum with respect to each dataset in the - plane, and contours are then drawn at , , and , corresponding to , , and . The minimum values of and the corresponding model parameters are given in table 2, which shows that the fit results are relatively insensitive to the mediator mass (the fits to the Fermi data display a mild preference for heavier mediators). Reasonably good fits to the Fermi and CCW data sets are obtained, with

(10) |

whereas the fit to the Daylan et al. data is poor. The data are compared to the simulated observed spectrum from the GC in fig. 2 for representative values of and , taking a mediator mass of .

The previous results are based upon the assumption of eq. (4) for the DM velocity disperion in the MW. The effect of using higher , using eq. 5, is shown in fig. 3, which results in somewhat lower central values of for the cross section and GeV for the mass.

CCW () | Fermi () | Daylan () | |||||||
---|---|---|---|---|---|---|---|---|---|

12 | 29.7 | 68 | -20.0 | 24.9 | 109 | -19.9 | 54.1 | 56 | -19.4 |

20 | 29.9 | 70 | -19.9 | 23.7 | 116 | -19.9 | 65.3 | 62 | -19.3 |

30 | 29.9 | 76 | -19.9 | 22.7 | 128 | -19.9 | 71.3 | 67 | -19.3 |

50 | 30.6 | 88 | -19.8 | 22.0 | 146 | -19.8 | 76.8 | 76 | -19.2 |

### ii.2 Limits from Dwarf Spheroidal Galaxies

An upper limit on the gamma-ray flux from DM annihilation in 18 dwarf spheroidal galaxies with kinematic data has been determined using Fermi-LAT data Ackermann:2013yva (). This can be used in conjunction with the -factors presented in table 1 to obtain an upper limit on . The strongest such constraint comes from the dwarf galaxy Draco. At a distance of 80 kpc and with a relatively large -factor and high velocity dispersion, it would be the most likely to exhibit signs of -wave annihilating dark matter.

Ackermann et al. give the combined limit on (annihilation into ) at 95% C.L. for 15 dwarf spheroidal galaxies. In our model, DM annihilates to , leading to a different gamma-ray spectrum, but in this section and the next we assume the resulting limit on the annihilation rate in both cases is approximately the same. (Note that the total energy deposition in the two cases is the same.)

The previously derived limit assumes -wave annihilation and therefore cannot be directly converted into a limit from -wave annihilation, as the different velocity dispersions of the dwarf spheroidals would have to be taken into account individually. If, however, we make the simplifying assumption that all dwarf spheroidal galaxies have velocity dispersions equal to the greatest value (that of Draco, with km/s), we can then use equation 3 to directly convert the limit to one on . This will lead to a constraint that is slightly more stringent than the true value, but sufficient for our purpose of showing that there is no tension with the GCE. The resulting upper limit on as a function of is shown in figure 4, along with the GCE best-fit regions. The weaker CMB constraint from energy injection at recombination Diamanti:2013bia () (also discussed in section IV.2) is also indicated there.

We see that the assumption of -wave annihilation rather than -wave completely eliminates the tension between the dwarf spheroidal constraints and the GCE. The former are softened by a factor of relative to the GCE signal. The constraints depend on the velocity dispersion profile assumed for the dwarfs, but even taking into account the uncertainties, the limiting cross section from dwarf spheroidal galaxies is far above the values required to explain the GCE.

### ii.3 Galaxy Cluster Limits

Searches for gamma rays from galaxy clusters can place more stringent constraints on our scenario. Although dwarf spheroidal constraints were weakened due to their smaller velocity dispersion, the converse is true for clusters: their larger velocity dispersions amplify the signal from -wave annihilating dark matter, relative to smaller systems.

Observations of the Coma Rephaeli:2015nca () and Virgo Ackermann:2015fdi () clusters have recently been analyzed by the Fermi-LAT Collaboration. The first of these references gave no limits on annihilating dark matter, while the second did so for -wave annihilations. We therefore derive the bound on -wave annihilating DM arising from the latter. For this purpose we adopt a value for the velocity dispersion of km/s for Virgo Girardi:1995iy ().

Limits on are derived at 95% C.L. for the Virgo cluster in Ackermann:2015fdi (), using a background model taking into account all Fermi 2-year catalog point sources as well as diffuse galactic and extragalactic spectra. We have converted them directly into limits on using equation 3, with one caveat: dark matter substructure—subhalos residing within the larger host halo—is expected to significantly boost the signal strength from -wave dark matter annihilation over what would be expected from the host halo alone. The constraints in ref. Ackermann:2015fdi () for the more conservative limit given assume a boost factor of from the substructure of the cluster. The substructure is not expected to have the same velocity dispersion as the host halo however, making the simple rescaling described in the previous section inapplicable for -wave annihilation. Subhalos generically have a significantly smaller velocity dispersion than the host halo, due to the fact that the velocity dispersion depends on the total mass of the subhalo where the dark matter is virially bound, not on that of the host halo. This can be seen in simulations such as RHAPSODY Wu:2012wu (); Wu:2012dg (), in which the number of galaxy cluster subhalos is found to drop off sharply with increasing maximum circular velocity (a power law index of -2.95) with no subhalos exceeding a third of the host halo’s maximum circular velocity. The contribution to the signal from subhalos is therefore weakened due to the velocity dependence of the annihilation cross section, offsetting the gains that come from the increased dark matter density. Ultimately, we choose a conservative approach and rescale the limits from Ackermann:2015fdi () by a factor of to remove the boost from the substructure for a self-consistent limit. The upper limits on are shown in fig. 4.

Similar limits have been found for several other clusters, including Coma and Fornax, using earlier Fermi data Ackermann:2010rg (). The Fornax cluster was subsequently reanalyzed with specific attention to the effects of subhalos and contraction due to baryonic infall Ando:2012vu (), leading to a more stringent upper bound on . As with the Virgo cluster, from this work we use the conservative limits neglecting the effect of substructure, which in ref. Ando:2012vu () are given alongside the more optimistic limits. We convert the constraints on Coma Ackermann:2010rg () (which does not account for substructure) and Fornax Ando:2012vu () directly into limits on , using velocity dispersions of km/s Girardi:1995iy () and km/s Drinkwater:2000dr () respectively; these are also included in fig. 4.

Although our best-fit parameters are consistent with older bounds from the Virgo and Fornax clusters, more recent observations of the Coma cluster are expected to give more stringent constraints due to its high dark matter density and larger velocity dispersion. Currently there are no limits on dark matter annihilation rates from the more recent observations, and such a study is beyond the scope of the present work.

### ii.4 Isotropic Gamma-Ray Background

The isotropic gamma-ray background (IGRB) could further constrain our scenario. As part of the DM annihilation contribution to this signal could be from even larger halos than the ones surrounding the galaxy clusters we considered in the previous section, it is possible that it could be enhanced if the annihilation cross section is velocity dependent. The most recent measurements of the IGRB can constrain the -wave annihilation cross-section to for conservative limits and for more optimistic limits corresponding to our adopted best-fit value of Liu:2016ngs (). Converting these limits into constraints on is not a simple matter, as arriving at an expected IGRB signal requires taking into account how the velocity dispersion varies for halos of different sizes and at different redshifts. Such a detailed analysis is beyond the scope of this work but would be interesting for future investigation.

## Iii Relic abundance from decaying dark matter

An immediate problem with -wave annihilating DM in the galactic center is that the corresponding cross section in the early universe would have been orders of magnitude greater, due to the larger relative velocities, leading to a highly suppressed relic density. The form of the Boltzmann equation which describes the time evolution of the number density for Dirac dark matter is

(11) |

where is the number density of a particle in thermal equilibrium with the photon bath. The equation for the evolution of the number density of the antiparticle is of the same form. We assume that there is no asymmetry between and , and therefore the total number density is given by

(12) |

Following the procedure of ref. Kolb:1990vq (), an approximate solution of the Boltzmann equation for the relic density is given by

(13) |

where , is the freeze-out temperature, and the effective degrees of freedom and are evaluated at . The thermally averaged cross section takes the form ; hence and for our -wave annihilation scenario where

(14) |

An approximate solution for is given by

(15) | |||||

Our fiducial fit, eq. (10), implies

(16) |

to be compared to the observed value Ade:2015xua (). Hence the thermally produced abundance is approximately 7000 times too small; we need a nonthermal production mechanism.

### iii.1 Decaying dark matter

A conceptually simple solution, similar to the superWIMP model proposed in Feng:2003xh (), is to suppose that today’s dark matter is the product of a heavier metastable state , that decayed into at temperatures below freeze-out of annihilations. For , this occurs at according to (16). Hence we need for to have a lifetime exceeding s. Such long lifetimes are suggestive of an analog of weak interactions in the dark sector. We consider representative effective interactions giving rise to decays , or , of the form

(17) |

where are heavy scales. Each operator is also accompanied by its Hermitian conjugate, which leads to decays of . These decay channels are chosen to illustrate constraints that can arise from big bang nucleosynthesis (BBN) and the cosmic microwave background (CMB). An alternate channel would be safe from these constraints. The decay rates corresponding to the first two operators are given by

(18) |

where the mass splitting is considered to be much less than , but greater than for decays into electrons. (We ignore phase space effects in the small region of parameter space where .) For the third operator, we are interested in larger mass splittings since must be at least . We use numerical results for its decay rate where needed. A fairly good fit is given by where for in GeV units, .

To obtain the relic density of the parent particle , we assume for definiteness an effective interaction

(19) |

giving rise to , where can be a light fermion of the standard model or in a hidden sector. The annihilation cross section for is

(20) |

To determine the relic density in this scenario, we again use eqs. 13 and III but now with , since the vector current operators of eq. (17) lead to -wave annihilation. The ultimate relic density of particles is related to the prior abundance of by . Curves of constant in the - plane for GeV are shown in figure 5. Here we consider two different scenarios: annihilations to electrons and positrons and to quark-antiquark pairs. Large mass splittings GeV lead to a reduction in that must be compensated by reducing the cross section by increasing . These estimates assume that coannihilations as well as inelastic scatters are unimportant for determining the DM relic density. This will be true (as we explore in detail in the following subsections) as long as , which is also consistent with the need for to be relatively long-lived. For small GeV, the desired relic density for and is independent of and requires GeV when and couple to and GeV when they couple to .

### iii.2 Coannihilations

Coannihilation processes can reduce the relic density of , which was assumed to be a small effect in the previous treatment. When the splitting between and is small, leading to , the effect can be estimated by replacing in eq. 12 with Griest:1990kh ()

(21) |

Here X represents any standard model particle, so the first term in the above equation is the total annihilation cross section and the second is the total coannihilation cross section. Eq. (12) with this effective cross section only describes the number density of until the freeze-out of this species, since after that point decays and inelastic scatterings can have a significant impact on the number density. In this low mass splitting limit, the relevant coannihilation cross section for the 4-fermion operator in eq. (17) has the same form as eq. (20), while the dipole operator gives Weiner:2012cb ()

(22) |

where is the electric charge of the fermions in the final state and is the fine-structure constant.

In the scenario where the relic density is determined entirely by coannihilation processes, i.e., when the operator in eq. (19) is not present, the correct relic density requires GeV or GeV. These are lower bounds, since increasing the strength of coannihilation processes would lead to underproduction of DM, while the larger relic density induced by decreased coannihilation can be offset by increased annihilation.

The resulting limits are shown in fig. 6. Decays of to -quarks require a relatively large mass splitting and consequently a more sophisticated calculation than the one described here, but the limits on from suppressing coannihilations are greatly subdominant to those from demanding that decays after freeze-out. We also note that the operators in eq. (17) lead to additional annihilation processes from the ones we have considered above, including for the four-fermion operator as well as and for the magnetic dipole operator. We have checked that these are negligible when the other constraints considered are satisfied.

### iii.3 Inelastic scattering

A further requirement for consistency of our relic density determination is that inelastic scatterings induced by the decay operators (17) are not important during the epoch between freeze-out and the significantly later freeze-out. Otherwise further depletion of the final abundance would occur due to scattering-induced transitions followed by annihilations. This leads to the criterion

(23) |

for the operator. We ignore the effect of transitions because the number density of relative to is extremely suppressed in our scenario at freeze-out. For GeV, freeze-out occurs at GeV, for which it is sufficient to compute the inelastic cross section in the elastic limit , and also approximating . We find that

(24) |

at the relevant energies. Performing the thermal average over electron energies gives , and we find from (23) the limit

(25) |

where is the abundance of 90 GeV DM relative to electrons at .

From the magnetic dipole operator, one has photon-mediated scattering from all charged particles that are in equilibrium at GeV, which we take to be plus their antiparticles. The cross section has a logarithmic infrared divergence in the limit from low-angle scattering. For , it is regulated more effectively by Debye screening than by the small value of , giving

(26) | |||||

For simplicity we cut off the divergence for all species using the Debye mass GeV. The thermal average of (26) is for the parameters of interest. The resulting bound analogous to (25) is

(27) |

No similar constraint arises for since quarks are not present in the plasma at temperature .

The bounds on and are shown in fig. 6. In both cases the limits derived from suppressing inelastic scattering are much stronger than those from suppressing coannihilation processes. This is because the number density of relativistic standard model scattering partners is much greater than the Boltzmann-suppressed number density of at freeze-out.

## Iv Constraints on decaying DM

To ensure that decays occur after freeze-out of -wave annihilations estimated in (16), we assume that for the relevant decay rate, with Hubble parameter and for GeV. Comparing to the decay rates (18), we obtain constraints on the parameter space in figs. 6,7, shown in the lower regions of the plots. In the unshaded central regions, decays occur after freeze-out and before big bang nucleosynthesis (BBN) or recombination. In the upper shaded regions, decays will disrupt BBN or the cosmic microwave background (CMB), due to the deposition of electromagnetic energy, as well as hadrons in the case of decays to , as we consider in the following subsections.

### iv.1 BBN constraints

For the first two operators of (17), leading to decays into photons or electrons, only the total energy deposited in the plasma is relevant for photoproduction or dissociation of light elements produced by BBN. We take the combined constraints from ref. Cyburt:2002uv () (see fig. 8 of that reference). An upper limit on as a function of lifetime is derived there, which we convert into a limit on as a function of , shown in fig. 6 for GeV. (The choice of determines .) Since the limit on is not monotonic in lifetime, BBN excludes a range of for a given value of .

The third operator of (17) leading to pairs
entails somewhat more
stringent constraints because of hadronic interactions that can
more efficiently disturb light element abundances Jedamzik:2006xz (). The limits depend not
only upon the total amount of energy deposited, but also the
energy per decay. By interpolating between the constraints of
Jedamzik:2006xz () calculated for different masses of decaying DM,
we find the BBN lower limit on versus shown
in fig. 7.^{5}^{5}5The relevant constraints are
inferred from figs. 9-10 of Jedamzik:2006xz (), in the region
s, which is insensitive to uncertainties in the
observed Li/Li abundance. The role of DM mass in that
reference (where the DM particle was assumed to decay completely to standard model particles)
is played by in the present context.

### iv.2 CMB constraints

For lifetimes s, electromagnetic energy deposition starts to distort the cosmic microwave background, superseding BBN constraints. We have computed the Planck-projected upper limits on the injected energy fraction as a function of lifetime using the tools of ref. Slatyer:2012yq () (see also ref. Cline:2013fm ()), where transfer functions are provided for computing the efficiency of energy deposition as a function of redshift for injections of photons or electrons at . For , can be used directly since the spectrum is monochromatic. For , must be convolved with the normalized energy spectrum of electrons from the 3-body decay, which in the limit of takes the form , where . Converting the limits on versus into the - plane results in the excluded regions shown in fig. 6. These extend to lifetimes greater than the age of the universe, not of interest in the present context, since would still be the principal component of the dark matter.

Projected Planck limits on the lifetimes for decays into have been given in ref. Cline:2013fm () for several DM masses. Interpolating those results we translate them into 95% C.L. limits on as a function of , shown in fig. 7.

### iv.3 Direct detection

For keV, it is possible to have direct detection through inelastic scattering on nuclei, . This is relevant for the magnetic dipole operator for which such small mass splittings are in the allowed region of fig. 6. We have roughly indicated the region excluded by direct searches there by taking the scattering rate to scale as where is the minimum velocity for an inelastic transition. It is given in terms of the DM-nucleus reduced mass as . Therefore the experimental limit on takes the form for some reference mass splitting, which we estimate to be keV by comparison to recent constraints on magnetic inelastic dark matter found in ref. Barello:2014uda (). The coefficient corresponds to the limit from elastic scattering (), which we take to be GeV by rescaling the constraints on dipolar dark matter from ref. Sigurdson:2004zp () according to the latest limits from the LUX experiment Akerib:2015rjg ().

In principle, the four-fermion operators in (17) could give rise to inelastic scattering on nucleons, by forming a loop from the electrons or quarks and considering virtual photon exchange between the loop and protons in the nucleus (see fig. 8). However the scattering rate is negligible since the required mass splitting or is too large to be excited in direct detection experiments. For smaller , there is an electron-loop mediated decay (decays into 1 and 2 photons are forbidden by gauge invariance or Furry’s theorem), but this is too slow to be of interest for , since the lifetime exceeds s and puts the model into the CMB-excluded region.

### iv.4 Fermi gamma ray line search

The magnetic dipole operator in eq. (17) induces annihilation to monochromatic gamma rays through -channel exchange of a particle. The cross section for this process has been calculated in Weiner:2012cb ():

(28) |

The Fermi-LAT collaboration has searched for such signals of DM annihilation in the Milky Way halo Ackermann:2015lka (). In the left plot in figure 6 we show the limits at 95% C.L. on the magnetic dipole operator from their search, assuming that the DM density follows a generalized NFW profile with and corresponding to a region of interest of around the galactic center to maximize the expected signal Ackermann:2015lka (). The line search limit is GeV, roughly equivalent to the bound we obtained from coannihilations.

### iv.5 Collider constraints

The operators we consider are also constrained by collider searches. For the operator in eq. (17), the relevant limits come from LEP, where the characteristic signature is missing energy and a photon which is radiated off the initial or . For our fiducial case of GeV, DELPHI monophoton searches constrain GeV at 90% C.L. Fox:2011fx (). For the magnetic dipole operator, the current most stringent constraint is from LHC monojet searches Barger:2012pf () requiring GeV at 95% C.L., a limit which is only slightly more constraining than searches for monophotons at LEP Fortin:2011hv () or the LHC Barger:2012pf (). The collider-disfavored region for the magnetic moment operator is more strongly excluded in our scenario by direct detection and the lifetime constraints. The operator is in principle limited by LHC monojet searches, but the small -quark content of the proton makes such limits very weak.

All of the exclusions discussed here were derived under the assumption that or collisions lead to stable final-state dark sector particles. Although it is possible for the produced in these collisions to decay inside the detector, in the regions of parameter space for which the collider limits are relevant, the mass splitting is so small that the softness of the decay products would render them undetectable.

For the operator of eq. (19), the relevant limits are again from LEP monophoton searches when and LHC monojet searches when . Limits from an ATLAS monojet search Aad:2015zva () as well as that from the previously mentioned DELPHI monophoton search are shown in fig. 5. In either case the correct relic density is compatible with current collider limits.

## V Mediator couplings

A large -wave cross section would generically run afoul of direct and indirect detection constraints if the DM coupled directly to SM particles. On the other hand, annihilation to light mediators , that subsequently decay into SM particles, can avoid this problem. If couples to as , the resulting -wave cross section at low velocities is given by

(29) |

An uncomfortably large coupling would be needed to match the fit to the GC excess.

However smaller values of can be sufficient if the cross section is Sommerfeld enhanced, which can naturally occur if the mediator is light. Defining , an analytic approximation to the enhancement factor is given by Cassel:2009wt ()

(30) |

for partial wave scattering, where

(31) |

with velocity in the center of mass frame. For a -wave process we take . is nonvanishing in the limit , so that velocity suppression of the -wave cross section is still present despite the enhancement, and has quasiperiodic resonant behavior as a function of .

The enhancement factor depends on the relative velocity of the particles, which in principle must be averaged over phase space. Ignoring the radial dependence of the annihilation cross-section, we can find an estimate of the average enhancement, which is given by

(32) |

using eq. (1). However for the parameter values of interest, we find that the dependence on is very weak and one can simply use - km/s with negligible error from omitting the average. To match the desired value of the cross section in (10), we need

(33) |

This relation as a function of is shown as the dashed line in fig. 9, while the analytic approximation for is indicated by the curves for two different values of the mediator mass.

By comparing to the required cross section, eq. (10), we find that the coupling constant can be reduced to a more comfortable value of . It can be somewhat further reduced by taking larger values of , as can be seen in figure 9; the left panel shows decreasing with , while the right shows increasing with . It was recently pointed out Blum:2016nrz () that approximations to the enhancement factor such as (30) may fail to satisfy partial wave unitarity in the resonant regions. We have checked that we are very far from any such violation however, for the parameter values of interest.

Finally, it has been noted in An:2016kie () that it is possible for two DM particles to capture into a bound state and then annihilate to mediators. The bound state formation process dominantly occurs in the -wave, so, if possible, it dominates over the direct -wave annihilation to mediators. In forming a bound state a mediator is emitted, so the mass of the mediator must be less than the binding energy of the ground state for this to occur, i.e.

(34) |

For the values of and that we consider to explain the GC excess, GeV for a bound state to form. In this work we have only considered mediator masses above this limit.

## Vi Conclusion

We have presented a scenario in which -wave annihilating dark matter could have significant indirect signals from the galactic center despite having a velocity-suppressed cross section. Although our immediate motivation was to reconcile a dark matter interpretation of the observed GC gamma ray excess with conflicting constraints from dwarf spheroidal galaxies, the framework presented here could be of more general interest.

Our key idea is to assume that the current generation of -wave annihilating DM is the decay product of a metastable predecessor particle , which has a thermal origin and decays after -wave annihilations of the stable DM have frozen out. This allows a large range of lifetimes that depend upon the - mass splitting and the mass scale of the effective interaction which leads to the decay. A number of constraints must be satisfied, including those coming from BBN, the CMB, direct detection, photon line searches, mono-X searches at colliders, and scattering in the early universe (which could deplete the DM abundance). They depend strongly upon the nature of the decays, which we have illustrated using the examples of , , and , but in all cases there is a significant region of the - parameter space in which all of the requirements can be satisfied.

The annihilation cross sections of interest for explaining the galactic center excess are larger than would generically occur in the presence of -wave suppression because of the low DM velocity in the GC. We nevertheless demonstrated a working example using reasonable coupling strengths, where the DM annihilates into light scalar mediators that mix with the Higgs boson and subsequently decay into quarks. We have shown that such models give a reasonably good fit to the observed GCE, while satisfying constraints from dwarf spheroidals by a comfortable margin.

Collider tests of our scenario are currently weaker than the consistency requirement that inelastic scatterings on standard model particles do not deplete the DM relic density in the early universe (due to strong -wave annihilations of ). For a narrow window of mass splittings MeV, direct SM searches provide a possible means of detection in the case of magnetic inelastic transitions.

Fermi observations of gamma rays from galaxy clusters may provide a more sensitive test of our scenario, due to the large velocity dispersion in clusters. We have shown that limits on DM annihilation from the Virgo cluster, while significantly stronger than limits from dwarf spheroidals, still are far from being in tension with this interpretation of the GCE. We hope that our work will motivate further studies of limits on DM annihilation in the Coma cluster, which has the potential to be more constraining because of its relatively high velocity dispersion.

Acknowledgments. We thank Kimmo Kainulainen, Maxim Pospelov, Subir Sarkar, Tracy Slatyer, Alex Drlica-Wagner and Wei Xue for helpful discussions or correspondence. We acknowledge support of the McGill Space Institute and the Natural Sciences and Engineering Research Council of Canada. We thank Niels Bohr International Academy for its generous hospitality during the initial stages of this work.

## References

- (1) L. Goodenough and D. Hooper, (2009), 0910.2998.
- (2) D. Hooper and L. Goodenough, Phys. Lett. B697, 412 (2011), 1010.2752.
- (3) D. Hooper and T. Linden, Phys. Rev. D84, 123005 (2011), 1110.0006.
- (4) K. N. Abazajian and M. Kaplinghat, Phys. Rev. D86, 083511 (2012), 1207.6047, [Erratum: Phys. Rev.D87,129902(2013)].
- (5) B. Zhou et al., Phys. Rev. D91, 123010 (2015), 1406.6948.
- (6) F. Calore, I. Cholis, and C. Weniger, JCAP 1503, 038 (2015), 1409.0042.
- (7) T. Daylan et al., Phys. Dark Univ. 12, 1 (2016), 1402.6703.
- (8) S. Murgia, presented at Fifth Fermi Symposium, 20- 24 Oct. 2014, http://fermi.gsfc.nasa.gov/science/mtgs/ symposia/2014/program/08 Murgia.pdf .
- (9) Fermi-LAT, M. Ajello et al., Astrophys. J. 819, 44 (2016), 1511.02938.
- (10) S. K. Lee, M. Lisanti, B. R. Safdi, T. R. Slatyer, and W. Xue, Phys. Rev. Lett. 116, 051103 (2016), 1506.05124.
- (11) R. Bartels, S. Krishnamurthy, and C. Weniger, Phys. Rev. Lett. 116, 051102 (2016), 1506.05104.
- (12) R. M. O’Leary, M. D. Kistler, M. Kerr, and J. Dexter, (2015), 1504.02477.
- (13) T. Linden, Phys. Rev. D93, 063003 (2016), 1509.02928.
- (14) I. Cholis, D. Hooper, and T. Linden, JCAP 1506, 043 (2015), 1407.5625.
- (15) E. Carlson and S. Profumo, Phys. Rev. D90, 023015 (2014), 1405.7685.
- (16) J. Petrović, P. D. Serpico, and G. Zaharijaš, JCAP 1410, 052 (2014), 1405.7928.
- (17) I. Cholis et al., JCAP 1512, 005 (2015), 1506.05119.
- (18) D. Gaggero, M. Taoso, A. Urbano, M. Valli, and P. Ullio, JCAP 1512, 056 (2015), 1507.06129.
- (19) W. de Boer, I. Gebauer, S. Kunz, and A. Neumann, (2015), 1509.05310.
- (20) E. Carlson, T. Linden, and S. Profumo, (2016), 1603.06584.
- (21) Q. Yuan and K. Ioka, Astrophys. J. 802, 124 (2015), 1411.4363.
- (22) S. K. Lee, M. Lisanti, and B. R. Safdi, JCAP 1505, 056 (2015), 1412.6099.
- (23) Fermi-LAT, M. Ackermann et al., Phys. Rev. Lett. 115, 231301 (2015), 1503.02641.
- (24) K. N. Abazajian and R. E. Keeley, Phys. Rev. D93, 083514 (2016), 1510.06424.
- (25) DES, K. Bechtol et al., Astrophys. J. 807, 50 (2015), 1503.02584.
- (26) DES, A. Drlica-Wagner et al., Astrophys. J. 813, 109 (2015), 1508.03622.
- (27) A. Geringer-Sameth et al., Phys. Rev. Lett. 115, 081101 (2015), 1503.02320.
- (28) D. Hooper and T. Linden, JCAP 1509, 016 (2015), 1503.06209.
- (29) S. Li et al., Phys. Rev. D93, 043518 (2016), 1511.09252.
- (30) DES, Fermi-LAT, A. Drlica-Wagner et al., Astrophys. J. 809, L4 (2015), 1503.02632.
- (31) M. R. Buckley et al., Phys. Rev. D91, 102001 (2015), 1502.01020.
- (32) B. Bertoni, D. Hooper, and T. Linden, JCAP 1512, 035 (2015), 1504.02087.
- (33) J. Liu, N. Weiner, and W. Xue, JHEP 08, 050 (2015), 1412.1485.
- (34) F. Calore, I. Cholis, C. McCabe, and C. Weniger, Phys. Rev. D91, 063003 (2015), 1411.4647.
- (35) M. Kaplinghat, T. Linden, and H.-B. Yu, Phys. Rev. Lett. 114, 211303 (2015), 1501.03507.
- (36) E. Hardy, R. Lasenby, and J. Unwin, JHEP 07, 049 (2014), 1402.4500.
- (37) Y. Zhao, X.-J. Bi, H.-Y. Jia, P.-F. Yin, and F.-R. Zhu, Phys. Rev. D93, 083513 (2016), 1601.02181.
- (38) M. Abdullah et al., Phys. Rev. D90, 035004 (2014), 1404.6528.
- (39) G. Elor, N. L. Rodd, and T. R. Slatyer, Phys. Rev. D91, 103531 (2015), 1503.01773.
- (40) B. Robertson and A. Zentner, Phys. Rev. D79, 083525 (2009), 0902.0362.
- (41) M. G. Walker et al., Astrophys. J. 667, L53 (2007), 0708.0010.
- (42) G. Battaglia et al., Mon. Not. Roy. Astron. Soc. 364, 433 (2005), astro-ph/0506102, [Erratum: Mon. Not. Roy. Astron. Soc.370,1055(2006)].
- (43) W. Dehnen, D. McLaughlin, and J. Sachania, Mon. Not. Roy. Astron. Soc. 369, 1688 (2006), astro-ph/0603825.
- (44) J. F. Navarro, C. S. Frenk, and S. D. M. White, Astrophys. J. 462, 563 (1996), astro-ph/9508025.
- (45) E. Romano-Diaz, I. Shlosman, Y. Hoffman, and C. Heller, Astrophys. J. 685, L105 (2008), 0808.0195.
- (46) F. Chen, J. M. Cline, A. Fradette, A. R. Frey, and C. Rabideau, Phys. Rev. D81, 043523 (2010), 0911.2222.
- (47) M. Cirelli and J. M. Cline, Phys. Rev. D82, 023503 (2010), 1005.1779.
- (48) Fermi-LAT, M. Ackermann et al., Phys. Rev. D89, 042001 (2014), 1310.0828.
- (49) M. I. Wilkinson et al., (2006), astro-ph/0602186, [EAS Publ. Ser.20,105(2006)].
- (50) M. I. Wilkinson et al., Astrophys. J. 611, L21 (2004), astro-ph/0406520.
- (51) M. G. Walker et al., Astron. J. 131, 2114 (2006), astro-ph/0511465, [Erratum: Astron. J.132,968(2006)].
- (52) A. Koch et al., Astrophys. J. 657, 241 (2007), astro-ph/0611372.
- (53) M. G. Coleman, K. Jordi, H.-W. Rix, E. K. Grebel, and A. Koch, Astron. J. 134, 1938 (2007), 0708.1853.
- (54) K. B. Westfall et al., Astron. J. 131, 375 (2006), astro-ph/0508091.
- (55) V. Belokurov et al., Astrophys. J. 647, L111 (2006), astro-ph/0604355.
- (56) J. D. Simon and M. Geha, Astrophys. J. 670, 313 (2007), 0706.0516.
- (57) M. Geha et al., Astrophys. J. 692, 1464 (2009), 0809.2781.
- (58) B. Willman et al., Astron. J. 142, 128 (2011), 1007.3499.
- (59) A. Cuoco, B. Eiteneuer, J. Heisig, and M. Krämer, (2016), 1603.08228.
- (60) C. Evoli, D. Gaggero, D. Grasso, and L. Maccione, JCAP 0810, 018 (2008), 0807.4730.
- (61) A. W. Strong, I. V. Moskalenko, and O. Reimer, Astrophys. J. 537, 763 (2000), astro-ph/9811296, [Erratum: Astrophys. J.541,1109(2000)].
- (62) G. Di Bernardo, C. Evoli, D. Gaggero, D. Grasso, and L. Maccione, JCAP 1303, 036 (2013), 1210.4546.
- (63) M. Cirelli et al., JCAP 1103, 051 (2011), 1012.4515, [Erratum: JCAP1210,E01(2012)].
- (64) P. Ciafaloni et al., JCAP 1103, 019 (2011), 1009.0224.
- (65) P. Agrawal, B. Batell, P. J. Fox, and R. Harnik, JCAP 1505, 011 (2015), 1411.2592.
- (66) J. M. Cline, G. Dupuis, Z. Liu, and W. Xue, Phys. Rev. D91, 115010 (2015), 1503.08213.
- (67) R. Diamanti, L. Lopez-Honorez, O. Mena, S. Palomares-Ruiz, and A. C. Vincent, JCAP 1402, 017 (2014), 1308.2578.
- (68) Fermi-LAT, M. Ackermann et al., Astrophys. J. 819, 149 (2016), 1507.08995.
- (69) Fermi-LAT, M. Ackermann et al., Astrophys. J. 812, 159 (2015), 1510.00004.
- (70) M. Girardi et al., Astrophys. J. 457, 61 (1996), astro-ph/9507031.
- (71) H.-Y. Wu, O. Hahn, R. H. Wechsler, Y.-Y. Mao, and P. S. Behroozi, Astrophys. J. 763, 70 (2013), 1209.3309.
- (72) H.-Y. Wu, O. Hahn, R. H. Wechsler, P. S. Behroozi, and Y.-Y. Mao, Astrophys. J. 767, 23 (2013), 1210.6358, [Astrophys. J.767,23(2013)].
- (73) M. Ackermann et al., JCAP 1005, 025 (2010), 1002.2239.
- (74) S. Ando and D. Nagai, JCAP 1207, 017 (2012), 1201.0753.
- (75) M. J. Drinkwater, M. D. Gregg, and M. Colless, Astrophys. J. 548, L139 (2001), astro-ph/0012415.
- (76) W. Liu, X.-J. Bi, S.-J. Lin, and P.-F. Yin, (2016), 1602.01012.
- (77) E. W. Kolb and M. S. Turner, Front. Phys. 69, 1 (1990).
- (78) Planck, P. A. R. Ade et al., (2015), 1502.01589.
- (79) J. L. Feng, A. Rajaraman, and F. Takayama, Phys. Rev. Lett. 91, 011302 (2003), hep-ph/0302215.
- (80) K. Griest and D. Seckel, Phys. Rev. D43, 3191 (1991).
- (81) N. Weiner and I. Yavin, Phys. Rev. D86, 075021 (2012), 1206.2910.
- (82) R. H. Cyburt, J. R. Ellis, B. D. Fields, and K. A. Olive, Phys. Rev. D67, 103521 (2003), astro-ph/0211258.
- (83) K. Jedamzik, Phys. Rev. D74, 103509 (2006), hep-ph/0604251.
- (84) T. R. Slatyer, Phys. Rev. D87, 123513 (2013), 1211.0283.
- (85) J. M. Cline and P. Scott, JCAP 1303, 044 (2013), 1301.5908, [Erratum: JCAP1305,E01(2013)].
- (86) G. Barello, S. Chang, and C. A. Newby, Phys. Rev. D90, 094027 (2014), 1409.0536.
- (87) K. Sigurdson, M. Doran, A. Kurylov, R. R. Caldwell, and M. Kamionkowski, Phys. Rev. D70, 083501 (2004), astro-ph/0406355, [Erratum: Phys. Rev.D73,089903(2006)].
- (88) LUX, D. S. Akerib et al., Phys. Rev. Lett. 116, 161301 (2016), 1512.03506.
- (89) Fermi-LAT, M. Ackermann et al., Phys. Rev. D91, 122002 (2015), 1506.00013.
- (90) P. J. Fox, R. Harnik, J. Kopp, and Y. Tsai, Phys. Rev. D84, 014028 (2011), 1103.0240.
- (91) V. Barger, W.-Y. Keung, D. Marfatia, and P.-Y. Tseng, Phys. Lett. B717, 219 (2012), 1206.0640.
- (92) J.-F. Fortin and T. M. P. Tait, Phys. Rev. D85, 063506 (2012), 1103.3289.
- (93) ATLAS, G. Aad et al., Eur. Phys. J. C75, 299 (2015), 1502.01518, [Erratum: Eur. Phys. J.C75,no.9,408(2015)].
- (94) S. Cassel, J. Phys. G37, 105009 (2010), 0903.5307.
- (95) K. Blum, R. Sato, and T. R. Slatyer, JCAP 1606, 021 (2016), 1603.01383.
- (96) H. An, M. B. Wise, and Y. Zhang, (2016), 1606.02305.