Model Independent analysis of MeV scale dark matter: I. Cosmological constraints

TL;DR

This paper develops a model-independent method to analyze MeV-scale dark matter's cosmological constraints, revealing that a large parameter space remains viable despite existing bounds.

Contribution

It introduces a new methodology to determine effective couplings of MeV-scale dark matter to hadrons, enabling cosmological constraint analysis independent of specific models.

Findings

Cosmological and astrophysical observations constrain some parameter space.

A significant portion of MeV-scale dark matter parameter space remains viable.

Methodology allows for model-independent analysis of light scalar dark matter.

Abstract

Recent results from several direct detection experiments have imposed severe constraints on the multi-GeV mass window for various dark matter (DM) models. However, many of these experiments are not sensitive to MeV scale DM as the corresponding recoil energies are, largely, lower than the detector thresholds. We reexamine the light scalar DM in a model-independent approach. In this first of a two-part work, we develop an appropriate methodology to determine the effective coupling of such a DM to hadrons, thereby allowing for the determination of the corresponding annihilation rates. We find that while the parameter space can be constrained using cosmological and astrophysical observations, a significantly large fraction is still viable. In the companion paper, we study the sensitivity of both direct detection experiments as well as colliders to such a DM.

Tables

Table 1. Table 1: List of various matrix elements and their relation to form factors.

$J^P=0^{-}$ : ${ℬ}_{1,2}=\pi ,K,\eta ,\cdots$		${ℬ}_1=\pi$,  ${ℬ}_2=\rho$
$⟨0\left|\overline{q}q\right|ℬℬ⟩$	$F_s$	$⟨0\left|\overline{q}{\gamma }^{\sigma }q\right|\pi \rho ⟩$	$F_v$  ${\epsilon }_{\mu \nu \rho \sigma }{p}_3^{\mu }{p}_4^{\nu }{ϵ}^{\rho }$
$⟨0\left|\overline{q}{\gamma }^{\sigma }q\right|ℬℬ⟩$	$F_v{q}^{\sigma }$
${ℬ}_{1,2}=\rho :(q^{\sigma }=p_3-p_4$ and $P^{\sigma }=p_3+p_4)$
$⟨0\left|\overline{q}q\right|\rho \rho ⟩$		$\left(F_s^1{m}_{\rho }{g}^{\mu \nu }+F_s^2{p}_3^{\nu }{p}_4^{\mu }\right){ϵ}_{\mu }\left(p_3\right){ϵ}_{\nu }\left(p_4\right)$
$⟨0\left|\overline{q}{\gamma }_{5}q\right|\rho \rho ⟩$		$F_p{ϵ}^{\mu \nu \rho \sigma }{p}_{3\rho }{p}_{4\sigma }{ϵ}_{\mu }\left(p_3\right){ϵ}_{\nu }\left(p_4\right)$
$⟨0\left|\overline{q}{\gamma }^{\sigma }q\right|\rho \rho ⟩$		$\left[F_v^1{q}^{\sigma }+F_v^2\left(g^{\mu \sigma }{P}^{\nu }-g^{\nu \sigma }{P}^{\mu }\right)+F_v^3{p}_3^{\nu }{p}_4^{\mu }\frac{q^{\sigma }}{m_{\rho }^2}\right]{ϵ}_{\mu }\left(p_3\right){ϵ}_{\nu }\left(p_4\right)$
$⟨0\left|\overline{q}{\gamma }^{\sigma }{\gamma }_{5}q\right|\rho \rho ⟩$		$F_a{ϵ}^{\mu \nu \rho \sigma }\left(q_{\rho }+P_{\rho }\right){ϵ}_{\mu }\left(p_3\right){ϵ}_{\nu }\left(p_4\right)$

Table 2. Table 2: Here, the value of g ρ = 5.01 subscript 𝑔 𝜌 5.01 g_{\rho}=5.01 is calculated from the decay constant of the ρ 𝜌 \rho -meson while the decay rates of the given processes are taken from PDG PhysRevD.98.030001 . Note that the couplings marked by † † \dagger are determined using vector meson dominance in the respective processes.

Coupling Constant	Value	Process
$g_{{\omega }_{782}}^{†}$	3.4$g_{\rho }$	${\omega }_{782}\to e^{-}{e}^{+}$
$g_{{\omega }_{782}\pi \rho }^{†}$	10.1${\mathrm{GeV}}^{-1}$	${\omega }_{782}\to {\pi }^{+}{\pi }^{-}{\pi }^{\circ }$
$g_{\varphi }^{†}$	2.7$g_{\rho }$	$\varphi \to e^{-}{e}^{+}$
$g_{\varphi \pi \rho }$	0.2${\mathrm{GeV}}^{-1}$	$\varphi \to \pi {\omega }_{782}$
$g_{{\omega }_{1420}}^{†}$	0.9$g_{\rho }$	${\omega }_{1420}\to e^{-}{e}^{+}$  assuming Br. frac. = ${10}^{-5}$
$g_{{\omega }_{1420}\pi \rho }$	3.8${\mathrm{GeV}}^{-1}$	${\omega }_{1420}\to \pi \rho$  assuming Br. frac. = 70%

Equations

$$\begin{array}[]{lcl}\mathcal{O}^{f}_{s}&=&\displaystyle\frac{\mathcal{C}^{f}_{s}}{\Lambda}\varphi^{\dagger}\varphi\;\bar{f}f\\[8.61108pt] \mathcal{O}^{f}_{p}&=&\displaystyle\frac{\mathcal{C}^{f}_{p}}{\Lambda}\varphi^{\dagger}\varphi\,\bar{f}\gamma^{5}f\\[8.61108pt] \mathcal{O}^{f}_{v}&=&\displaystyle\frac{\mathcal{C}^{f}_{v}}{\Lambda^{2}}i~{}(\varphi^{\dagger}\partial_{\mu}\varphi-\partial_{\mu}\varphi^{\dagger}\varphi)\;\bar{f}\gamma^{\mu}f\\[8.61108pt] \mathcal{O}^{f}_{a}&=&\displaystyle\frac{\mathcal{C}^{f}_{a}}{\Lambda^{2}}i~{}(\varphi^{\dagger}\partial_{\mu}\varphi-\partial_{\mu}\varphi^{\dagger}\varphi)\;\bar{f}\gamma^{\mu}\gamma^{5}f\\[8.61108pt] \mathcal{O}_{\gamma}&=&\displaystyle\frac{\mathcal{C}_{\gamma}}{\Lambda^{2}}(\varphi^{\dagger}\varphi)F_{\mu\nu}F^{\mu\nu}\\[8.61108pt] \mathcal{O}_{\tilde{\gamma}}&=&\displaystyle\frac{\mathcal{C}_{\tilde{\gamma}}}{\Lambda^{2}}(\varphi^{\dagger}\varphi)F_{\mu\nu}\tilde{F}^{\mu\nu}\ ,\end{array}$$

$$\frac{dn}{dt}+3H(t)n=-⟨\sigma v⟩(n^2-n_{eq}^2),$$

$$g_{\rho }=\sum _{\mathrm{bosons}}{g}_{\mathrm{bosons}}{\left(\frac{T_{\mathrm{bosons}}}{T_{\gamma }}\right)}^4+\frac{7}{8}\sum _{\mathrm{fermions}}{g}_{\mathrm{fermions}}{\left(\frac{T_{\mathrm{fermions}}}{T_{\gamma }}\right)}^4.$$

$$g_s=\sum _{\mathrm{bosons}}{g}_{\mathrm{bosons}}{\left(\frac{T_{\mathrm{bosons}}}{T_{\gamma }}\right)}^3+\frac{7}{8}\sum _{\mathrm{fermions}}{g}_{\mathrm{fermions}}{\left(\frac{T_{\mathrm{fermions}}}{T_{\gamma }}\right)}^3.$$

$$\frac{dY}{dx}=\frac{m_{\phi }^3⟨\sigma v⟩}{H(m_{\phi })x^2}\left(Y_{eq}^2-Y^2\right),$$

$$H(T=m_{\phi })=\sqrt{\frac{4{\pi }^{3}G{g}_{\rho }}{45}}{m}_{\phi }^2.$$

$$Y_{eq}=\frac{45}{2{\pi }^4}\sqrt{\frac{\pi }{8}}\frac{g_{\phi }}{g_{\rho }}{x}^{3/2}{e}^{-x}$$

$$Y_{\infty }\simeq x_{f}H(m_{\phi })/m_{\phi }^3⟨\sigma v⟩.$$

$${\rho }_{\phi }=m_{\phi }{Y}_{\infty }{T}_0^3{\left(\frac{a_f{T}_f}{a_0{T}_0}\right)}^3.$$

$${\Omega }_{\phi }{h}^2=\sqrt{\frac{4{\pi }^{3}G{g}_{\rho }(x_f)}{45}}\frac{x_f{T}_0^3{g}_0}{{\rho }_c⟨\sigma v⟩g_{\rho }(x_f)}.$$

$$\mathcal{M}=⟨{\mathcal{B}}_1{\mathcal{B}}_2|{\mathcal{O}}_I|\phi \phi ⟩=\frac{{\mathcal{C}}_I^f}{\Lambda }⟨{\mathcal{B}}_1{\mathcal{B}}_2|J_I^{\bar{f}f}|0⟩⟨0|J_{\phi }^I|\phi \phi ⟩.$$

$${\mathcal{O}}_{eff.}=\sum _{i=c,b}\frac{{\alpha }_s}{12\pi \Lambda m_i}{\phi }^{†}\phi G^{a\mu \nu }{G}_{\mu \nu }^a,$$

$${\theta }_{\mu }^{\mu }=-\frac{9{\alpha }_s}{8\pi }{G}^{a\mu \nu }{G}_{\mu \nu }^a+\sum _{\text{light}}{m}_q\bar{q}q,$$

$${\mathcal{O}}_{eff.}=\sum _{i=c,b}\frac{2}{27\Lambda m_i}{\phi }^{†}\phi \left[-{\theta }_{\mu }^{\mu }+\sum _{\text{light}=u,d,s}{m}_q\bar{q}q\right].$$

$${\mathcal{O}}_s=\frac{{\phi }^{†}\phi }{\Lambda }\left(\sum _{i=c,b}\frac{2}{27m_i}\left[-{\theta }_{\mu }^{\mu }+\sum _{q=u,d,s}{m}_q\bar{q}q\right]+\sum _{q=u,d,s}\bar{q}q\right).$$

$$F_s=-\sum _{k=c,b}\frac{2}{27m_k}{\theta }_{\mu }^{\mu }+\left(\frac{2}{27m_c}+\frac{2}{27m_b}\right){\Gamma }_m+\Gamma +\left(\frac{2}{27m_c}+\frac{2}{27m_b}\right){\Delta }_m+\Delta$$

$$\begin{array}[]{rcl}\Gamma_{m}&=&\displaystyle\langle{\cal B\,B}|m_{u}\bar{u}u+m_{d}\bar{d}d|0\rangle\\[4.30554pt] \Gamma&=&\displaystyle\langle{\cal B\,B}|\bar{u}u+\bar{d}d|0\rangle\\[4.30554pt] \Delta_{m}&=&\displaystyle\langle{\cal B\,B}|m_{s}\bar{s}s|0\rangle\\[4.30554pt] \Delta&=&\displaystyle\langle{\cal B\,B}|\bar{s}s|0\rangle\\[4.30554pt] \theta&=&\displaystyle\langle{\cal B\,B}|\theta_{\mu}^{\mu}|0\rangle\ .\end{array}$$

$$\begin{array}[]{rcl}\displaystyle q^{\mu}\langle 0|\bar{q}\gamma_{\mu}{q}|{\cal B\,B}\rangle&=&\displaystyle F_{v}\,q^{2}\\[4.30554pt] \displaystyle\langle 0|m_{q}\,\bar{q}{q}|{\cal B\,B}\rangle&=&\displaystyle F_{v}\,q^{2}\\[4.30554pt] \displaystyle\langle 0|\bar{q}\gamma_{\mu}{q}|{\cal B\,B}\rangle&=&\displaystyle\frac{(F_{s}^{\text{light}}+F_{s}^{\text{heavy}})q_{\mu}}{q^{2}}\ ,\end{array}$$

$$F_v=\frac{1}{q^2}\left(\frac{2}{9}\theta +\frac{7}{9}{\Gamma }_m+\frac{7}{9}\Delta \right).$$

$$\mathcal{L}=\frac{f_{\pi }^2}{4}\text{Tr}{\partial }_{\mu }U{\partial }^{\mu }{U}^{†}+\frac{B{f}_{\pi }^2}{2}\text{Tr}(M^{†}U+U^{†}M).$$

$$\pi(x)\equiv\frac{1}{\sqrt{2}}\,\left(\begin{array}[]{cc}\pi^{0}&\sqrt{2}\pi^{+}\\ \sqrt{2}\pi^{-}&-\pi^{0}\end{array}\right)\,$$

$$\pi(x)\longrightarrow\left(\begin{array}[]{ccc}\displaystyle\frac{\pi^{0}}{\sqrt{2}}+\frac{\eta}{\sqrt{6}}&\pi^{+}&K^{+}\\[6.45831pt] \pi^{-}&\displaystyle\frac{-\pi^{0}}{\sqrt{2}}+\frac{\eta}{\sqrt{6}}&K^{0}\\[6.45831pt] K^{-}&\displaystyle\overline{K^{0}}&\displaystyle\frac{-2\,\eta}{\sqrt{6}}\end{array}\right)\ .$$

$$\bar{q}q=-\frac{\partial {\mathcal{L}}_{\mathrm{QCD}}}{\partial m_q}=\frac{\partial {\mathcal{H}}_{\mathrm{QCD}}}{\partial m_q}.$$

$$\begin{array}[]{rcl}\displaystyle\frac{f_{\pi}^{2}}{2}\text{Tr}(M^{\dagger}U+U^{\dagger}M)&=&\displaystyle(m_{u}+m_{d})\pi^{+}\pi^{-}+\frac{(m_{u}+m_{d})}{2}\pi_{0}^{2}\\[8.61108pt] &+&\displaystyle(m_{u}+m_{s})K^{+}K^{-}+(m_{d}+m_{s})K^{0}\bar{K^{0}}+\frac{(m_{u}+m_{d}+4m_{s})}{6}\eta^{2}\\[8.61108pt] &+&\displaystyle\frac{(m_{u}-m_{d})}{\sqrt{3}}\pi^{0}\eta.\end{array}$$

$$\begin{array}[]{rcl c rcl}\displaystyle m_{\pi}^{2}&=&\displaystyle B(m_{u}+m_{d}),&&\displaystyle m_{K^{+}}^{2}&=&\displaystyle B(m_{u}+m_{s})\\[4.30554pt] \displaystyle\quad m_{K^{0}}^{2}&=&\displaystyle B(m_{d}+m_{s})&&m_{\eta}^{2}&=&\displaystyle B\frac{(m_{u}+m_{d}+4m_{s})}{3}.\end{array}$$

$$\displaystyle\Big{\langle}\pi\pi\Big{|}\sum_{q=u,d}\bar{q}q\Big{|}0\Big{\rangle}=B.$$

$$m_{\pi }^2=\hat{m}B,\hat{m}\equiv m_u+m_d,$$

$$\displaystyle\Big{\langle}\pi\pi\Big{|}\sum_{q=u,d}\bar{q}q\Big{|}0\Big{\rangle}=\frac{m_{\pi}^{2}}{\hat{m}}\ .$$

$$\begin{array}[]{rcl c rcl}\displaystyle\Big{\langle}K^{+}K^{-}\Big{|}\sum_{q=u,d}\bar{q}q\Big{|}0\Big{\rangle}&=&\displaystyle\frac{m_{\pi}^{2}}{2\,\hat{m}}\ ,&&\displaystyle\langle K^{+}K^{-}|\bar{s}s|0\rangle&=&\displaystyle\frac{2m_{K^{+}}^{2}-m_{\pi}^{2}}{2\,m_{s}}\ ,\\[6.45831pt] \displaystyle\Big{\langle}\eta\eta\Big{|}\sum_{q=u,d}\bar{q}q\Big{|}0\Big{\rangle}&=&\displaystyle\frac{m_{\pi}^{2}}{3\,\hat{m}}\ ,&&\displaystyle\langle\eta\eta|\bar{s}s|0\rangle&=&\displaystyle\frac{3\,m_{\eta}^{2}-m_{\pi}^{2}}{4\,m_{s}}\ ,\end{array}$$

$${\theta }_{\mu }^{\mu }=-{\partial }_{\mu }\pi {\partial }^{\mu }\pi +2m_{\pi }^2\pi \cdot \pi$$

Full text

Keywords: Dark Matter, MeV scale scalar, Belle-II

1 Introduction

The last half a century has witnessed the accumulation of overwhelming evidence for gravitational interactions between visible (and stable) particles and non-luminous matter on a multitude of scales, from the galactic to the cosmological. Starting with rotation curves in spiral galaxies Rubin:1970zza , gravitational lensing measurements Moustakas:2002iz ; 2012A&A…545A..71V , recent observations of cluster collisions (Bullet Cluster) Clowe:2006eq , a temperature anisotropy in the spectrum of Cosmic Microwave Background Radiation Hinshaw:2012aka ; Abadi:2002tt ; Ade:2015xua , there are a large variety of observations for which the Dark Matter (DM) hypothesis provides the most compelling explanation.

However, all these observations are indirect and, thus, it is difficult to ascertain whether it is a particulate DM that underlies these anomalies or whether the latter are but manifestations of our lack of understanding of gravity at different cosmological scales. These issues have, rightly, been explored and several alternates to standard gravity proposed Milgrom:2001ny ; Famaey:2011kh . On the other hand, it should be noted that even if such modifications of gravity do exist, theories incorporating these will not be necessarily unrelated from a model involving particles as DM Calmet:2017voc .

Despite the lack of direct evidence (i.e., one achieved under controlled conditions), it is the notion of particle DM that has seen the most development. The reasons are twofold. For one, the experimental discovery of particle DM is more viable in comparison to the verification of modified gravity models. Moreover, DM candidates arise naturally in a host of particle physics models that are motivated primarily to address other issues that are unanswered by the Standard Model (SM). To substantiate such models, a variety of experiments have been proposed, with many being already under operation (or even having outlived their use). In principle, this could be done in three ways: a) Satellite based indirect detection experiments like Fermi-LAT Ackermann:2015zua , PAMELA Adriani:2008zr , AMS Lin:2015taa , depend on the annihilation of a pair of DM particles into SM particles which can produce rare antimatter cosmic rays (positrons, anti-protons or antideuterons), neutrinos, monochromatic photons or continuous $\gamma$-ray spectrum. Although there, occasionally, have been claims of anomalies in the data, unfortunately the experiments have failed to validate each other’s putative positive sightings, resulting in further constraints. b) Direct detection experiments that, typically, identify the nuclear recoils produced by the scattering between DM and the detector’s (target) nuclei. c) Collider searches based on the production of DM look for the excesses (over the SM expectations) in final states with large missing momentum. (Note, though, that collider experiments are indicative at best, for these can only verify the stability of the putative DM candidate over detector dimensions, not over cosmological timescales.)

Taking cue from recent null results in the LUX Akerib:2016vxi , PandaX-II Tan:2016zwf and XENON100 Aprile:2016swn experiments for $m_{DM}>$ 6 GeV, we concentrate, here, on MeV scale DM(m_{DM}\mathrel{\hbox to0.0pt{\raise 2.15277pt\hbox{<}\hss}{\lower 2.15277pt\hbox{\sim}}}3\,{\rm GeV}). Light DM can easily evade many of these direct and indirect detection experiments because of the low momentum transfer (lower than the threshold)222However, if the DM family has light as well as heavy DM, then it can be detected in (direct) detections as in the case of boosted DM. In this case, DM can be energetic if heavy DM decays into the lighter one. Such DM candidates are not included in our analysis, as the details of the particular model are paramount in such scenarios, whereas we attempt only a model-independent analysis., and, consequently, the lack of signals in this range has motivated discussions of a DM with a particularly low mass Alexander:2016aln . Similarly, to explain perceived anomalies in the 511 keV $\gamma$-rays observed by INTEGRAL, the cosmic $\gamma$-ray background at 1-20 MeV and the details of large scale structure, quite a few models Boehm:2003bt ; Ahn:2005ck ; Borodatchenkova:2005ct ; Hooper:2007tu with a light DM were invoked. Again, WIMPless DM also accommodates DM masses in the MeV scale Feng:2008ya . These models emerge naturally from gauge-mediated supersymmetry breaking where DM naturally satisfies the current relic density without its mass and interaction being restricted to the weak scale.

All of the above mentioned models are well motivated but lack experimental support. Moreover, in case the DM particle is the only new particle (in the dark sector) within the reach of a particular experiment while other new species are much heavier, it will be very difficult to distinguish the underlying theories. Therefore, model-independent studies of the DM are powerful as they proffer a way ahead without being constrained to a very specific scenario.

The purpose of this article is to study the parameter space for MeV scale DM in a model-independent way. We first construct effective operators describing interactions between a scalar DM particle and the (visible) SM sector. The constraints on these, as obtained from DM relic abundance, CMB and the counting of relativistic degrees of freedom are then analysed.

2 Higher Dimension operators

The interaction of the dark matter with the SM sector is completely unknown barring, of course, the gravitational one. All that we know is that the DM does not interact either strongly or electromagnetically333While models have been proposed wherein the DM does have a very tiny charge, these tend to be baroque, and do not fit within well-motivated scenarios going beyond the SMOverduin:2004sz .. On the other hand, if the DM does not interact at all with the SM particles (except gravitationally), then there would, essentially, be no way to directly confirm their existence. More importantly, with the DM particles having been produced profusely during the post-inflation reheating phase and shortly thereafter, without such interactions, the relic density today would tend to be too large, thereby more than over-closing the universe, an eventuality that can be avoided only by tuning the initial conditions.

The interaction that, thus, must be posited could be in the form of a detailed and ultraviolet-complete model (such as that in the minimum supersymmetric standard model) or in the shape of an effective field theory. It is the latter approach that we adopt here, choosing to profess an ignorance of the underlying theory (the UV-completion). In other words, we would augment the SM with the DM particle and posit that the latter interacts with the known particles through certain higher-dimensional operators (without any explicit mediator being considered). With the typical energy scale of the processes under consideration being much smaller than the dominant mass scale of the theory, such an approach is irreproachable.

Our assumption, thus, is that the only new relevant field is the scalar444The results are identical for a pseudoscalar DM., with all other new species being too heavy to be relevant in the contexts of both terrestrial experiments/observations as well as the cosmological evolution of the relic density. Since we are interested in a DM with a mass of at most a few GeVs, the only relevant SM states are the photon and the gluon, the leptons (including neutrinos) and the quarks of the first two generations. The bottom-quark is, at best, only marginally relevant. Consequently, we consider operators including this limited set of particles alone. Furthermore, to be consistent with low-energy constraints, we do not admit flavour changing operators.

Assuming $SU(3)\otimes U(1)_{\rm em}$ symmetry555All but ${\cal O}^{f}_{s,p}$ respect the full SM gauge symmetry. These two too can be altered trivially to respect the full symmetry, albeit at the cost of introducing an extra factor of $\langle H\rangle/\Lambda$, where $H$ is the SM Higgs doublet. For the present analysis, this distinction is essentially irrelevant., the operators for a complex scalar field666It should be noted here that analogous results can be achieved for a real scalar as well. are

[TABLE]

where $f$ is an arbitrary SM fermion and $\Lambda$ is the scale of new physics. Note that we could also write operators akin to $\mathcal{O}_{\gamma}$ and $\mathcal{O}_{\tilde{\gamma}}$, but for the gluons instead. As for the ${\cal C}$’s (the dimensionless Wilson coefficients corresponding to the various operators), we would be normalizing these to either zero or unity (denoting the absence or presence of the said operator). The results will, thus, depend on the mass of DM and the the scale $\Lambda$. Indeed, with each operator presumably arising from a specific DM-SM interaction in a UV-complete theory, the conclusions reached from an analysis such as ours can be easily rescaled to obtain constraints on the parameters of the underlying theory.

3 Relic Abundance

The model independent framework developed in the preceding section allows us to constrain the parameter space for any MeV-scale spin-0 Dark Matter candidate. Particular attention needs to be paid to the constraints from the relic abundance, the cosmic microwave background and, on account of the lightness, that from the counting of relativistic degrees of freedom. In this section, we consider only the first, relegating the others to a later section.

3.1 The formalism

We restrict our discussions to the context of DM that had, primarily, been produced thermally and was in equilibrium with the SM sector. The relic abundance of non-thermal DM, on the other hand, depends crucially on the conditions when it was produced. This, being intricately tied to the specifics of the dynamics can only be addressed within the context of a particular model, and, hence, does not fall under the ambit of a model-independent analysis such as ours. Since the WIMPs are presumed to be produced thermally, the relic abundance calculationSteigman:2012nb ; Yu:2011by can proceed as usual while taking care of some subtleties owing to the small mass. For the WIMP to stay in thermal equilibrium, it needs to interact with the SM sector, with the strength(s) being sufficiently large to beat the expansion rate of the universe.

During its evolution, the (stable) spin-0 particle $\varphi$ (of mass $m_{\varphi}$) was in thermal equilibrium until a certain epoch. Similarly, the SM particles (barring, possibly, the neutrinos) were also in thermal equilibrium with the photon gas. The latter determines the temperature of the thermal soup, and this we shall denote by $T_{\gamma}$. The evolution of $\varphi$ is given by the Boltzmann equation, namely

[TABLE]

where $n$ is the number density of $\varphi$ ($n_{eq}$ being its equilibrium value), $H$ is the Hubble expansion rate and $\langle\sigma v\rangle$ is the thermally averaged cross-section for DM annihilation.

Before attempting a general solution of eq.(2), let us consider some general properties. For any massive particle, the number density at equilibrium depends on the ratio of its mass and the temperature of the plasma, $x\equiv m_{\varphi}/T$. For a stable particle that is relativistic yet at equilibrium ($x\ll 1$), the annihilation processes as well as pair production are proceeding at comparable rates. The consequent equilibrium density is given by $n=3\zeta(3)g_{\varphi}T^{3}/(4\pi^{2})$ where $g_{\varphi}$ denotes its degrees of freedom. On the other hand, if it is nonrelativistic, i.e., its mass is much larger than the ambient energy ($x\gg 1$), the plasma does not have sufficient energy to drive pair production; yet, the pair annihilation proceeds. Consequently, its equilibrium abundance falls exponentially as the temperature drops below the mass of the particle, yielding $n=g_{\varphi}(m_{\varphi}T/(2\pi))^{3/2}e^{(-m_{\varphi}/T)}$.

Applying the above to the WIMP, in the immediate aftermath of its production, it would have been in equilibrium due to the balance between its interactions with the SM particles as its interactions with this sector were strong enough to beat the expansion rate. If it continues to be in equilibrium, then as the universe cools to $T\ll m_{\varphi}$, the WIMP would become nonrelativistic and its abundance today would have been negligible. However the very structure of its interactions (essentially the higher-dimensional nature of the couplings), stipulates that it must fall out of equilibrium for $T\ll\Lambda$. Naively, it might seem that this would have occurred when the WIMP was still ultrarelativistic, since $m_{\varphi}\ll\Lambda$. On the other hand, N-body simulationsAbadi:2002tt for structure formation requires the DM to be non-relativistic777In this scenario, DM perturbations grew in the matter dominated era forming a gravitational well. (This could not have been initiated by ordinary matter as it could not have clustered due to the radiation pressure.) Ordinary matter could now fall into this well, thereby allowing an early start of the structure formation and formation of the fine structure in universe.. In other words, $\varphi$ should decouple from the thermal soup only when it had become non-relativistic in the radiation dominated era. This could happen if, thanks to the exponential suppression of $n_{eq}$, the WIMPs became so rare that the interaction rate fell below the expansion rate. No longer affected by interactions, these fall out of equilibrium with the abundance freezing out, (i.e, their number in a comoving volume becomes constant).

The freeze-out temperature $T_{f}$, namely that at the epoch when the DM number density freezes out, can be determined in terms of the mass and the interaction strengths, as we discuss now. In the radiation-dominated regime, it is useful to express the total energy density in terms of the photon energy density thereby defining $g_{\rho}$, the effective number of degrees of freedom associated with total energy density, namely

[TABLE]

In addition to this, the entropy($S$) in a comoving volume ($S=sa^{3}$) is a conserved quantity which enables us to define the number density in terms of the “yield” $Y=n/s$ and also, analogously, introduce $g_{s}$, the number of effective degrees of freedom associated with the entropy:

[TABLE]

With the density of a nonrelativistic species falling faster, $g_{\rho}$ and $g_{s}$ differ noticeably only when there are relativistic particles present that are not in equilibrium with photons. Within the SM, this occurs for neutrinos. Similarly, for m_{\varphi}\mathrel{\hbox to0.0pt{\raise 2.15277pt\hbox{<}\hss}{\lower 2.15277pt\hbox{\sim}}}6 MeV, the DM will contribute to relativistic degrees of freedom and, thence, entropy. For such light DM, the explanation of the relic abundance by way of the DM having reached thermal equilibrium with the SM sector is excluded by current observations which we discuss in the next subsection. Although there are models which use a different line of approach, namely asymmetric or non-thermal Iminniyaz:2011yp ; Dev:2013yza , we shall desist from doing so, and will no longer consider this range. Therefore, we can safely assume $g_{\rho}\simeq g_{s}$ in our analysis.

Entropy conservation also implies $a(t)T=$ constant (here, $a(t)$ is the scale factor of the universe) and, hence, $dT/dt=-H(t)T$. Effecting a change of variables, $T\to x\equiv m_{\varphi}/T$, we, then, have the famous Boltzmann equation, viz.

[TABLE]

where $H(m_{\varphi})$ is the Hubble expansion rate (in the radiation dominated universe) calculated at the epoch when the temperature equals $m_{\varphi}$ and is given by888In the radiation dominated universe, the scale factor $a(t)$ goes as $t^{1/2}$, while the temperature-time relation is given by $t^{2}\,T^{4}=45/8\pi^{3}G$. Together, this gives $H(T)$.

[TABLE]

A caveat needs to be entered here. A key ingredient in reaching eq.(3) is the assumption that entropy is conserved throughout the era of interest. However, this statement may not necessarily be true for MeV-range dark matter, especially if $500\,{\rm MeV}\leq m_{\varphi}\leq 1\,{\rm GeV}$ (see Fig. 7). The freeze out temperature (equivalently, $x_{f}\equiv m_{\varphi}/T_{f}$) for this range may lie around the QCD phase transition999See Sec.6.1 for a derivation of $x_{f}$.. Consequently, the entropy may not be conserved at this epoch101010It has been suggested Aoki:2006we , though, that the QCD transition in the early universe is not a real phase transition but an analytic cross-over. As we argue next, the distinction is of no consequence in the current context.. However, as the entropy is overwhelmingly determined by the contributions from relativistic particles such as $e^{\pm}$, $\gamma$ and $\nu$’s, this small possible non-conservation can be neglected altogether.

To obtain the present density of DM particles, we need to solve eq.(3) in terms of the final freeze out abundance $Y_{\infty}$ (at $x=\infty$). While this, unfortunately, can be done only numerically, it is instructive to consider an approximate analytic solution. Before freeze-out, $Y_{\varphi}$ was close to its equilibrium value, $Y_{eq}$

[TABLE]

which is exponentially suppressed111111Although, after freeze out, the abundance is larger than what its equilibrium value would have been, this approximation of $Y\approx Y_{eq}^{\rm freeze-out}$ is an excellent one, especially for understanding the structure of the solution.. Integrating eq.(3) from the freeze out temperature $x_{f}$ until very late times ($x=\infty$), we get

[TABLE]

With the energy density for the now non-relativistic DM, being given by $\rho_{\varphi}=m_{\varphi}\,n_{\varphi}$, post freeze-out, it simply falls as $a^{-3}$. Denoting the freeze-out epoch (i.e., when $Y$ has reached the asymptotic value $Y_{\infty}$) by the temperature $T_{f}$ and the scale factor $a_{f}$, with the corresponding quantities today being given by $T_{0}$ and $a_{0}$ respectively, we have $n(a_{f},T_{f})=Y_{\infty}\,T^{3}_{f}$, and, today,

[TABLE]

Simultaneously, the number of effective degrees of freedom changes from $g_{\rho}(x_{f})$ at the freeze out epoch to $g_{0}=3.36$ operative today, and $g_{0}\,a_{0}\,T^{3}_{0}=g_{\rho}(x_{f})\,a_{f}\,T^{3}_{f}$. It is customary to parametrize $\rho_{\varphi}\equiv\Omega_{\varphi}h^{2}\,\rho_{c}$, where $\rho_{c}=1.05375\times 10^{-5}\,h^{2}\left({\rm GeV}/c^{2}\right){\rm cm}^{-3}$ is the critical density of the universe with the Hubble constant today being expressed as $H_{0}=h\times 100$ km s*-1* Mpc*-1*. We have, then,

[TABLE]

Both the WMAP 2013ApJS..208…19H and the Planck Ade:2015xua satellite observations determine the relic density very well, with the latter suggesting $\Omega_{DM}h^{2}=0.1199\pm 0.0022$. Clearly, we must have $\Omega_{\varphi}\leq\Omega_{DM}$. This can be translated to constraints on the parameter space available to the theory. As the interaction between the DM and the SM sector increases, so does $\langle\sigma\,v\rangle$, and, as a consequence, the relic abundance of DM today decreases (see eq. 4).

3.2 Bounds from Relic abundance

The dimension-5 operators in eq.(1) would lead to $\sigma\propto\Lambda^{-2}$. On the contrary, for nonrelativistic particle, the dimension-6 operators stipulate $\sigma\propto m_{\varphi}^{2}\,\Lambda^{-4}$. Consequently, for the second set, a given value of $\Omega_{\varphi}h^{2}$ would require $m_{\varphi}$ to increase nearly quadratically with $\Lambda$. On the other hand, for the first set, the “right value” of $m_{\varphi}$ should have only a weak dependence on $\Lambda$.

To delineate the exact parameter space (as opposed to basing our conclusions on analytic results obtained from an approximations as we have been making so far), we implement these additional operators in micrOMEGAs 4.1Belanger:2014vza using FeynRulesAlloul:2013bka . To best understand the consequences, we only incorporate a single operator structure (from eq. 1) at a time. For the operators involving a fermionic current, we consider two different cases, namely $(a)$ one wherein all the SM fermions participate equally121212Clearly, given the mass range of the scalar, the third-generation quarks play little or no role., i.e., all the ${\cal C}^{f}$s of a given class are unity and $(b)$ the leptophilic case, namely where only the leptons participate (and equally) while the quarks do not. In Fig.1, we depict the contours in the $m_{\varphi}$–$\Lambda$ plane corresponding to $\Omega_{\varphi}h^{2}=0.1199$. The width in the contour due to the measurement error in $\Omega_{\varphi}h^{2}$ is virtually unobservable. In each case, the area below the curves would correspond to a larger annihilation cross-section (thanks to a smaller $\Lambda$) and, hence, a DM relic density smaller than what the Planck collaboration measures. In other words, this is the parameter space that is observationally allowed (with the remainder ostensibly being contributed by some other source). Note that, for the dimension-6 operators (Fig.1), the relation between $m_{\varphi}$ and $\Lambda$ is nearly quadratic, as expected. On the other hand, for the dimension-5 operators (Fig.1), $\Lambda$ increases much slower with $m_{\varphi}$. Understandably, the dimension-5 operators are sensitive to much larger values of $\Lambda$.

Expectedly, for the leptophilic DM, the value of $\Lambda$ for a given $m_{\varphi}$ is much lower than the case of the democratic coupling131313Here, we assume the simplistic point of view that the light quarks can be treated as (pseudo-)asymptotic states. This, of course, is untenable, and has been assumed only for illustrative purposes. We return to this point later.. This is but a reflection of the fact that, now, fewer final states are available to the DM annihilation process. And, finally, the dependence on the chirality structure is rather minimal, a reflection of the fact that the fermion masses in the problem are quite small.

The corresponding constraints for a real scalar can be divined from those discussed above by realizing that a complex scalar can be expressed as $\varphi_{c}=(\varphi_{1}+i\,\varphi_{2})/\sqrt{2}$. In other words, the relic density of the complex field may be expressed in terms of those for the real fields as $\Omega=\Omega_{1}+\Omega_{2}$. For identically parametrized effective Lagrangians, the vertices for the real scalar theory would have an extra factor of 2 (as compared to those for the complex field). On the other hand, one must account for identical fields in the final state when a pair of DM particles is being produced as a result of SM-pair annihilation. The consequent change in $T_{f}$ is, however, only a minor one (at the level of 10% or lower) and has little bearing on the relic abundance calculation. Now, $\Omega\propto\langle\sigma\,v\rangle^{-1}\sim\Lambda^{n}$, where $n=2\,(4)$ for dimension-5 (dimension-6) operators. The aforementioned factor of 2 in the vertex, along with a factor of half in the thermal average (owing to identical particles in the initial state), thus, implies that the constraint on $\Lambda_{{\rm real}\varphi}$ would be, approximately, a factor of $2^{3/2}\,(2^{3/4})$ stronger than those derived above for dimension-5 (dimension-6) operators. This is borne out by explicit calculations.

3.3 Caveat to the calculation

Until now, we have been considering the case where a pair of DM particles, annihilate into a pair of SM particles, treating the latter as asymptotic states. In other words, it was assumed, naively, that a DM-pair annihilating to hadrons could be well-approximated by $\varphi\,\varphi^{*}\rightarrow q_{i}\,\bar{q_{i}}$ with the quarks hadronizing subsequently. However, when the mass of the DM particle is of the order of a quark mass, the relative momentum between the quark-antiquark pair is small and a bound state ensues. This, obviously, would need a different calculational scheme. For the mass range that we are considering, this is of relevance only in regards to the three light quarks. In particular, note that such a DM species with interaction strengths that we are investigating would freeze out only around the QCD phase transition temperature (see Fig.7). This brings with the added complication that, even if the DM is considerably heavier than a light quark, the annihilation products would hadronize immediately (on the scale of the annihilation time), perhaps into a pair of bound states, thereby utterly disallowing the approximation of quarks as quasi-asymptotic states. Thus, if we want to admit unsuppressed DM-quark interactions, it is imperative that we consider annihilation to bound states, and we set up the formulation next.

4 DM annihilation to various bound states

For $m_{\varphi}>2m_{\pi}$ , the DM can annihilate to pions, kaons and other light mesons, with more and more channels being available to heavier DM particles. Nonetheless, if the coupling to $u/d$–quarks are unsuppressed, the dominant effect arises from DM annihilation to pions. The calculation of the relevant rates involves the determination of the matrix elements of the operators for hadronic states, and we begin by discussing these.

4.1 Matrix Element and Form factors

For arbitrary bound states ${\cal B}_{1,2}$, the matrix element for the process $\varphi+\varphi^{\dagger}\rightarrow{\cal B}_{1}+{\cal B}_{2}$, driven by an operator $\mathcal{O}_{I}$ listed in eq.(1), is given by

[TABLE]

where the operator has been factorized into a product of two currents. Since the DM, by definition, does not suffer strong interactions, the vacuum saturation approximation is almost exact, a result that we use in the second step above. Clearly, the valence quark content of the hadrons ${\cal B}_{1,2}$ must be conjugate of each other for the corresponding matrix element to be nonzero. In other words, together, they should be a flavour-singlet pair. For simplicity, now onwards, we will assume that the DM interacts universally with quarks and leptons, and set $C_{I}^{f}=1,\,\forall f$ (though only for a given current structure). Wherever it is pertinent, we shall indicate the difference that the relaxation of this assumption will entail.

The matrix elements can be parametrized in terms of form factors multiplying momentum-dependent structures dictated by Lorentz covariance, parity transformation and other symmetries, like isospin, wherever applicable. Some of these are listed in Table 1. Others can be parametrized analogously. Note that many final state are missing in Table 1. For example, the $\pi\,K$ final state is precluded (at least to the lowest order in electroweak theory) by the assumption that the DM interactions do not admit flavour violation. Similarly, our assumption of identical couplings of the DM to the up– and down-quarks implies isospin symmetry and, consequently, final states such as $\pi\eta$ are strongly suppressed141414If the couplings ${\cal C}_{s}^{u}$ and ${\cal C}_{s}^{d}$ are unequal, this suppression is inoperative and the $\pi\eta$ channel opens up, with an amplitude proportional to $({\cal C}_{s}^{u}-{\cal C}_{s}^{d})$. We, however, do not consider such an explicit $SU(2)$-violation any further.. Indeed, as we shall see shortly, for scalar currents, the amplitude for this final state is proportional to the quark mass difference $m_{u}-m_{d}$. For vector currents, on the other hand, the process suffers an additional $v^{4}$ suppression, as is expected for a pair of scalars annihilating to another through a vector current.

4.2 Scalar Form Factors

We begin by attempting to relate the simplest of the quark currents, viz. the scalar $\bar{q}\,q$ to mesonic currents. Naively, for the DM masses of interest here, couplings to heavy quarks should not play a role. However, they actually do, courtesy quantum corrections. For example, integrating out the heavy quarks would result in an effective operator of the form $\varphi^{\dagger}\varphi G^{\mu\nu}G_{\mu\nu}$. At the one-loop order, this can be estimated by calculating a triangle diagram with heavy quarks as propagators yielding

[TABLE]

where, for brevity’s sake, higher powers in $m_{i}^{-1}$ have been neglected. We have also explicitly omitted the top-quark contribution as it is highly suppressed at this scale.

The operator in eqn.(6), though, suffers higher order corrections, and an accurate perturbative calculation thereof is rather cumbersome. It is useful, however, to recast it in terms of the trace anomaly and appeal to the known renormalization group flow of the energy momentum tensor $\theta^{\mu\nu}$ Voloshin:1985tc ; Chivukula:1989ze . In the present context, the trace of the QCD $\theta^{\mu\nu}$ is given by

[TABLE]

and the operator in eqn.(6) becomes

[TABLE]

Finally, for the scalar operator, we have

[TABLE]

Note that, if the coefficients ${\cal C}_{s}^{q}$ were different from unity, the terms in the equation above would be trivially modified.

Using the result above, we can define the scalar form factor as

[TABLE]

where

[TABLE]

4.3 Vector Form Factors

For a real (i.e., one that carries no charge or any other additive quantum number) (pseudo)-scalar meson ${\cal B}$, a matrix element of the form $\langle 0|\bar{q}\gamma_{\mu}q|{\cal BB}\rangle$ would, necessarily, vanish identically. This would be the case for the $\pi^{0},\eta$, but not necessarily for the $K^{0}$. For the latter and for other charged mesons, the scalar and pseudo-scalar form factor can be related to vector and axial-vector form factor respectively through

[TABLE]

where $q^{\mu}=(p_{\bar{q}}-p_{q})^{\mu}$. Using expressions analogous to those in the preceding subsection to express contributions due to the heavy quarks, the total vector form factor and hence, matrix element for vector interaction is given by

[TABLE]

In other words, the vector form factor can be written in terms of the scalar form factors, and the extraction of the latter suffices. Therefore, in the rest of the paper, we will focus on ways to extract the scalar form factors.

An accurate determination of the form factors requires the expression of the quark-current in terms of hadronic currents and several approaches are possible. A particularly simple and elegant formalism is afforded by chiral perturbation theory ($\chi$PT) and the use of dispersion relations. Analogous techniques have been used in studying the decay of a light scalar into hadrons Donoghue:1990xh ; Truong:1989my ; Chivukula:1989ze ,and, in the next section, we adapt these to our case.

5 Using Chiral Perturbation Theory ($\chi$PT)

We begin by recapitulating the key results, derived within chiral perturbation theory ($\chi$PT), that are useful to our work. For an in-depth discussion of the subject, numerous resources  Leutwyler:1993iq ; Bijnens:1994qh ; Scherer:2002tk exist.

As is well known, $\chi$PT describes the low energy dynamics of QCD. In its simplest version, corresponding to the existence of just two massless quarks ($u,d$), QCD admits an exact $SU(2)\otimes SU(2)$ chiral symmetry, and the corresponding $\chi$PT lagrangian is described by

[TABLE]

Here, $f_{\pi}$ is the pion decay constant, and the dynamical degrees of freedom are encoded in a $2\times 2$ matrix $U\in SU(2)$, which can be parametrized as $U=e^{i\pi(x)/f_{\pi}}$, where

[TABLE]

transforming in the adjoint representation of $SU(2)$ would be identified with the physical pions. In eqn.(13), $B$ is an arbitrary coupling constant (whose physical significance is yet to be ascertained) while $M$ is a constant $2\times 2$ matrix to be related to masses. The aforementioned Lagrangian would exhibit $SU(2)_{L}\otimes SU(2)_{R}$ symmetry only if the matrix $M$ also transforms appropriately, viz., under the action of the chiral symmetry group, $U\rightarrow V_{L}\,U\,V^{\dagger}_{R}$ and $M\rightarrow V_{L}\,M\,V^{\dagger}_{R}$, where $V_{L,R}$ denote the respective transformations under the two $SU(2)$s. Given that $M$, unlike $U$, is not a dynamical variable, such a transformation may seem strange. However, note that, for massive quarks, QCD lagrangian does not admit axial symmetry. Indeed, even in the absence of $M$, the ground state of the above lagrangian is not symmetric under axial symmetry. In other words, the symmetry is spontaneously broken leading to three goldstone bosons with odd parity. Similarly, for unequal quark masses, the $SU(2)_{V}$ symmetry is also lost. Keeping this in mind, a nondynamic $M$ can be perceived as a perturbation that explicitly breaks $SU(2)_{A}$, thereby rendering the pions to be only pseudo-Goldstone bosons, as also breaking $SU(2)_{V}$ by a small amount proportional to the difference $m_{d}-m_{u}$. Consequently, matrix elements can be expanded in powers of the mass term, or, equivalently, as $\mathcal{O}(p^{2})$ corrections.

For three flavours, the symmetry is enlarged to $SU(3)_{L}\otimes SU(3)_{R}$, or, equivalently, to $SU(3)_{V}\otimes SU(3)_{A}$. Now $U\in SU(3)$ and $\pi(x)$ refers to the full pseudoscalar octet, viz.,

[TABLE]

Similarly, $M$, too, gets promoted to a $3\times 3$ matrix. This is the theory that we shall work with.

We now consider the quark operators of interest. As has been demonstrated in the preceding section, the vector form factors can be expressed in terms of the scalar ones. The operators $\bar{q}_{i}q_{i}$ (where $q_{i}$ denote the light quarks) can be expressed as

[TABLE]

On the other hand, the term, $\text{Tr}(M^{\dagger}U+U^{\dagger}M)$, in the $\chi$PT Lagrangian can be expanded up to the second order to obtain the pion, kaon and $\eta$ mass terms, namely

[TABLE]

Hence, the masses are given as

[TABLE]

Using the fact that the expectation values of the respective Hamiltonians for the two theory should be equal, we have

[TABLE]

On the other hand,

[TABLE]

and, therefore,

[TABLE]

In a similar vein, we have (neglecting the difference $m_{d}-m_{u}$)

[TABLE]

where the matrix elements for the $\bar{s}s$ current are obtained by differentiating ${\cal L}_{\rm QCD}$ with respect to $m_{s}$.

At this order, the trace of the energy momentum tensor reads

[TABLE]

Therefore, the corresponding form factors (see eqn.(10) for definition) are rendered

[TABLE]

and, consequently, the total matrix element becomes

[TABLE]

It should be noted that, in the limit of isospin symmetry ($m_{u}=m_{d}=\hat{m}/2$), the various scalar form factors obey

[TABLE]

Similar expressions hold for the kaonic form factors. Using these form factors and the expressions for the cross sections, we may obtain the relic abundance.

It should be realized, though, that the form factors derived so far have been defined within the lowest order $\chi$PT. This is reflected by the form factors being constants rather than functions of momenta. In other words, the aforementioned values only reflect the values of the form-factors at a particular momentum scale, defined by the decay/interaction which these are extracted from. As we shall shortly see, the higher-order corrections can be quite important. Consequently, we postpone the calculation of the relic abundance until after at least some of these corrections are evaluated.

The $\chi$PT Lagrangian can be expressed as a power series in the exchange momentum $p$; terms containing quark masses or external scalar or pseudo scalar fields are ${\cal O}(p^{2})$ whereas external vector or axial-vector fields are ${\cal O}(p)$. The NLO terms in the Lagrangian contains terms that are ${\cal O}(p^{4})$ or, in other words, suppressed by further factors of ${\cal O}(p^{2}/\Lambda_{\rm QCD}^{2})$. With $\Lambda_{\rm QCD}\sim 200\,{\rm MeV}$, clearly, a perturbative calculation of the higher-order effects is valid only for small momentum exchanges. In the present context, this translates to a limit on the dark matter mass, viz. m_{\varphi}\mathrel{\hbox to0.0pt{\raise 2.15277pt\hbox{<}\hss}{\lower 2.15277pt\hbox{\sim}}}300\,{\rm MeV} Chivukula:1989ze , for a perturbative expansion to make sense. Instead, we calculate the form factors using dispersion relations as this method relies solely on general principles and data.

5.1 Form factors from Dispersion Relations

We now discuss how the form factors can be extracted from scattering data using dispersion relations. This would allow us to determine the deviations, as compared to the preceding section, wrought in the relic abundances. Since the pion and kaon final states result in the dominant contributions, we would be concentrating primarily on these two form factors. In this, we largely follow the methodology developed in Refs.Donoghue:1990xh ; Ananthanarayan:2004xy .

5.1.1 From $\chi$PT to Dispersion Relations

With interactions amongst the hadrons switched on, at the one-loop level, diagrams as in Fig.2 would also contribute. While we have denoted only a subset of the one-loop diagrams, multiple intermediate states do contribute. And, with the hadron-hadron interactions being strong, there is no a priori compelling reason to limit ourselves to only one-loop results. In other words, to write down the S-matrix for such a system, we need to include contributions from channels like $\pi\pi\rightarrow\pi\pi$, $\pi\pi\rightarrow K^{+}\,K^{-}$, $\pi\pi\rightarrow 4\pi$, $\pi\pi\rightarrow\eta\,\eta$ etc. The direct calculation of the loops is, of course, a very difficult task. Instead, we take recourse to determining the imaginary part (wherever applicable) using the Cutkosky rules and, subsequently, calculating the real part using dispersion relations. Although, all channels do contribute to the total amplitude, in this work, we are restricting to two channels alone, viz. $\pi\pi$ and $K\bar{K}$ as these are expected to overwhelmingly dominate upto $m_{\varphi\varphi}=1.4\,{\rm GeV}$Moussallam:1999aq i.e., for dark matter of mass $\leq$ 700  MeV. Above $m_{\varphi\varphi}=1.4\,{\rm GeV}$, states such as $4\pi$ and $\eta\eta$ come into play, with the contribution of $f_{0}(1500)$, in particular, expected to be felt beyond $m_{\varphi\varphi}=1.5\,{\rm GeV}$ Moussallam:1999aq . Nonetheless, the 2-channel approximation is expected to be a very good one, as the other channels ($X$) mentioned above typically suffer from either kinematic restrictions (leading to a vanishing imaginary part) or small $\varphi\varphi\to X$ amplitudes151515In addition, for certain choices of $X$, even the amplitude for $\pi\pi$ to $X$ is suppressed as well.. We shall see this explicitly in results below.

Given that the initial state ($\varphi\varphi$) is well-described by $J^{P}=0^{+}$, we are interested only in final states with $I=0$ and $J^{P}=0^{+}$. Under the assumption that the other channels can be entirely neglected, the S-matrix for meson scatterings (which lives in the aforementioned subspace) can be further reduced to a $2\times 2$ unitary submatrix given by

[TABLE]

where $j,k=\pi\,,K$ and $M_{jk}$ is the corresponding element of the transition matrix.

Using the unitarity of the S-matrix, the imaginary part of the transition matrix can be expressed as

[TABLE]

and, similarly, for the imaginary part of form factor viz.,

[TABLE]

with $\Theta(x)$ being the well-known step function.

More explicitly, the Cutkosky rules determine the discontinuity, and, hence, the imaginary part of the scattering amplitude through the Schwarz reflection principle. In doing this, it needs to be realized thatDonoghue:1990xh the form factors of eqn.(10), treated as functions of $s\equiv p^{2}$, where $p_{\mu}$ is the momentum transfer, are analytic in the complex $s$-plane, except for a cut along the positive real axis. For the two-meson (${\cal BB}$) case, the cut starts at $s=4m_{\cal B}^{2}$. While for $s>4m_{\cal B}^{2}$, the identifications in eqn.(10) hold, for $0<s<4m_{\cal B}^{2}$, the function $\theta_{\cal B}(s)$ represents the matrix element $\langle 0|\theta_{\mu}^{\mu}|{\cal BB}\rangle$. Similarly, for real and negative values of $s$, it corresponds to $\langle{\cal B}|\theta_{\mu}^{\mu}|{\cal B}\rangle$ and is real. This, in turn, implies that the values above and below the cut are complex conjugates of one another.

Consider, for example, the region $4m_{\varphi}^{2}<s<4m_{K}^{2}$, the only allowed process for DM annihilation is, thus161616Here, we include the possibility that $\varphi$ may represent a real scalar. Similarly, $\pi\pi$ includes both $\pi^{+}\pi^{-}$ and $\pi^{0}\pi^{0}$. $\varphi\varphi^{*}\to\pi\pi$, with $3\pi$ final states being precluded by isospin invariance. Consequently, the only allowed final state rescattering is $\pi\pi\to\pi\pi$, and, in the limit of identical masses for the charged and neutral pions, is an entirely elastic process. Concentrating on the rescattering, the in- and out-states only differ in phase. For such a single-channel case, if we denote the form-factor on the upper side of the cut by $F^{\pi}_{s}$, then

[TABLE]

where $\delta_{\pi}(s)$ is the $I=0,\,J=0$ pion-scattering phase shift. Consequently, as the cut is approached from the above, $F^{\pi}_{s}\,\exp\left(i\,\pi\,\delta_{\pi}\right)$ is a real quantity.

Once $\delta_{\pi}$ is known (from data), what remains is to determine $F^{\pi}_{s}$. If, for $s\rightarrow\infty$, the phase shift $\delta_{\pi}$ tends to a finite value and $F^{\pi}_{s}$ does not grow faster than a power of $s$, then the form factor is known to be given by the Omnès functionOmnes:1958hv $\Omega(s)$ as

[TABLE]

where $P(s)$ is a polynomial to be fixed by the behaviour of $F_{s}(s)$. It is straightforward to prove that, for $\delta(s\rightarrow\infty)\rightarrow\alpha\pi$, the Omnès function is monomially suppressed, namely, $\Omega(s\rightarrow\infty)\rightarrow s^{-\alpha}$. Moreover, the high energy behaviour of QCD dictates that, asymptotically, the form factor behaves as $s^{-1}$. So, for pion-pion scattering, where phase shift at $s=\infty$ goes as $\pi$, the Omnès function would go as $1/s$ and $P(s)$, immediately, is constrained to be a constant, to be determined from the value of form factor at $s=0$.

For larger values of $s$, further channels come into play. Restricting ourselves, as argued for above, to two channels, the expression above is generalised to

[TABLE]

where the interaction amplitude for a process can be expanded in partial waves with angular momentum $l$:

[TABLE]

For scalar form factors, we need to consider only the $l=0$ part. Taking a cue from the single-channel case, the S-matrix for two channel process can be parametrized as

[TABLE]

where $\cos\theta$ determines the mixing between the two channels. Clearly, for $\cos\theta=1$, the two channels decouple entirely, and the solutions are to be obtained independently from $\pi\pi\to\pi\pi$ and $KK\to KK$ respectively, namely

[TABLE]

However, as is expected, and as can be ascertained from data (see, for example, Fig.4 of Ref.Au:1986vs ), $\cos\theta\neq 1$. Non-trivial values of $\cos\theta$ essentially parametrize the relative strength of, say, $KK$ admixture in the determination of the $\pi\pi\to\pi\pi$ scattering. Known as the elasticity parameter, $\cos\theta$ is a function of energy, angular momentum and isospin. Along with the phase shifts, we treat $\cos\theta$ as an experimental input.

It is useful to parametrize the consequent form factors through

[TABLE]

where the Clebsch-Gordan coefficient occurring in the projection of the $\pi\pi$ state with $I=0$ is shifted to $F^{K}$, and can be thought of as a relative normalisation.

5.1.2 The two channel solution: iterative procedure

The parameters of the S-matrix (or, equivalently, the T-matrix) may be determined using various scattering data and certain theoretical constraints in the Roy-Steiner equationsRoy:1971tc . (An easier route would be use the existing determination of phase shifts and inelasticity parameters, such as those found in Refs.Hoferichter:2012wf ; Au:1986vs . Before we describe the calculation of $\Omega^{i}_{j}$’s, let us express the form factors in terms of their values at low momentum transfers $(s=0)$, viz.

[TABLE]

Similarly, $\Omega^{i}_{j}$ are normalised as $\Omega^{1}_{1}(0)=\Omega^{2}_{2}(0)=1$ and $\Omega^{1}_{2}(0)=\Omega^{2}_{1}(0)=0$. The parameters $p_{i}$ are related to the slopes ($d\theta_{i}/ds$ and $d\Omega^{i}_{j}/ds$) evaluated at $s=0$. Relatable to $SU(3)$ breaking effects, and estimated by the use of unsubtracted dispersion relations, these lie in the range $(0.9,1.1)$. In other words, $p_{i}=1$ represents a very good approximation and reproduces the zeroth-order relations of eqn.(24). Armed with these “initial conditions”, we follow an iterative procedure similar to that outlined in Refs.Donoghue:1990xh ; Ananthanarayan:2004xy , to calculate the form factors. It is convenient to begin with, say, $\Delta_{\pi}$ and $\Delta_{K}$, as these can be handled in a fashion similar to the one channel case. We confirmed that, asymptotically these two solutions indeed go as $s^{-1}$. Furthermore, in the zeroth approximation, it is assumed that form factors behave as $F^{\pi}(s)=1$, $F^{K}(s)=\lambda$ where $\lambda$ is a real number. Then, the imaginary part of $F^{i}(s)$ at every iteration is computed via eqn. 29 and the real part is computed using

[TABLE]

The nested integral equations can be evaluated using standard numerical methods (see, e.g., Ref.Moussallam:1999aq ). The resultant form factors, as displayed in Fig. 3 and Fig. 3, are consistent with those obtained in Ref.Monin:2018lee ; Winkler:2018qyg .

5.2 Revisiting the Relic abundance for scalar interaction

Having obtained the form factors for scalar interaction, we have, once again, plotted the contours in the $m_{\varphi}-\Lambda$ plane in Fig.4 using form factors determined at the lowest order in $\chi$PT. Whereas Fig.1 corresponds to the evidently untenable assumption that the quarks may be treated as free particles, Fig.4 reflects the inclusion of bound state effects. One particular effect was expected. Contrary to the previous case, the annihilation channel into quarks open up only for $m_{\varphi}\sim m_{\pi}$. This was only to be expected as the pion is lightest of the hadrons. That values of $m_{\varphi}$ slightly smaller than $m_{\pi}$ are allowed too, is ascribed to the fact that the DM does have a non-zero momentum, albeit, somewhat smaller than its mass.

More interesting is the height of the peak. This sudden change is to be understood in terms of the structure of the annihilation cross sections for the different channels. For annihilations into a pair of leptons, the rate is proportional to $\Lambda^{-2}$. While this scaling originates from the structure of the effective Lagrangian in eqn.(1) and remains operative for annihilations into mesons as well, now one has to include the effect of the nontrivial wavefunction of the bound state. For a two-pion final-state, this appears in the matrix element as $F_{s}(q^{2})$ which has mass dimension 2, and the corresponding cross section would scale as $\sigma\propto m_{\pi}^{4}/(\Lambda^{2}\hat{m}^{2}\hat{s})$. With $\hat{s}\sim 4m_{\pi}^{2}$. it is the smallness of $\hat{m}$ that pushes up the cross section near the threshold. A similar bump exists for $m_{\varphi}$ slightly larger than $m_{K}$, but with a much smaller amplitude owing to the analogue of $\hat{m}$ being dominated by the much larger $m_{s}$. Even smaller is the contribution from the $\eta$-meson (as is testified to by Fig.4). However, note that, for $m_{\varphi}\sim 1\,{\rm GeV}$ the curve for the lowest-order bound-state analysis (wherein we have used only the $\pi\pi,KK,\eta\eta$ final states) is already veering very close to the free-quark one (which has been provided in Fig.4 for reference). Clearly the inclusion of more bound states will, asymptotically, render the contours very close to each other. For m_{\varphi}\mathrel{\hbox to0.0pt{\raise 2.15277pt\hbox{<}\hss}{\lower 2.15277pt\hbox{\sim}}}1\,{\rm GeV}, on the other hand, the true annihilation cross-sections are, typically, larger than what what transpires for free quarks, and, hence, the preferred value of $\Lambda$ somewhat larger.

We now move on, from the lowest order analysis, to the inclusion of higher orders, through the use of dispersion relations. Fig.5 displays the corresponding contours. As is evident from a comparison with Fig.4, the gross behaviour is quite similar. Overall, the increase in the sizes of the form factors, that the higher-order terms entail, results in a further rise in the preferred value of $\Lambda$. The bump around $m_{\varphi}\approx m_{K}$ is to be understood in terms of the interference between the pion and kaon form factors whose details are secured in the data of pion and kaon phase shifts.

5.3 Vector Interactions

For a final state comprised of a pair of (pseudo)scalars, the form factors for a vector current can be calculated in terms of those for the corresponding scalar current. As discussed in Sec.4.3, such annihilation cross sections suffer a $v^{4}$ suppression and, as a result, for vector interactions, such final states contribute very little (as compared to, say, the leptons). On the other hand, were the final state to comprise, say a pion and a rho meson or a pair of vector mesons, the cross sections would no longer be suppressed, and the conclusions would change drastically. In view of this, we next make an estimate of the time like form factors for these states.

In doing this, we would be using results pertaining to the well-studied electromagnetic current, which, for light mesons, reads

[TABLE]

Assuming isospin symmetry, we further have

[TABLE]

and, finally,

[TABLE]

where $\lambda$ denotes the polarization state of the $\rho$-meson and $F_{\rm em}(q^{2})$ is the electromagnetic form factor. The other form factors discussed in Table 1 can be defined analogously.

While, for the $\rho\rho$ final state, the results of Ref.deMelo:2016lwr can be used in a relatively straightforward manner, the $\pi\rho$ state needs more work, especially in determining the time-like electromagnetic form factor in the region q^{2}\mathrel{\hbox to0.0pt{\raise 2.15277pt\hbox{<}\hss}{\lower 2.15277pt\hbox{\sim}}}4\,{\rm GeV}^{2}. Towards this end, we make use of the vector meson dominance model,wherein, for q^{2}\mathrel{\hbox to0.0pt{\raise 2.15277pt\hbox{<}\hss}{\lower 2.15277pt\hbox{\sim}}}4\,{\rm GeV}^{2}, the major contributions accrue from the $\omega(782)$, $\phi(1020)$ and $\omega(1420)$.

Restricting ourselves, for the time being, to a single vector meson $\hat{V}_{\mu}$, the part of the Lagrangian governing its propagation, and the interactions with an external current (such as that corresponding to a pionic current) $J^{\mu}$ can be parametrized by Sakurai ; OConnell:1995nse

[TABLE]

Here, $\hat{V}_{\mu\nu}$ and $\hat{F}_{\mu\nu}$ are the usual field strengths corresponding to $\hat{V}_{\mu}$ and the photon field $A_{\mu}$. While $g_{{}_{VJ}}$ determines the strength of the coupling, the term containing $g_{{}_{V}}$ parametrizes the kinetic mixing (between the two vector fields) that is allowed in an abelian theory. The presence of this last term (which could also be thought of as a momentum-dependent $\gamma$–$\rho$ coupling vanishing at $q^{2}=0$) calls for a re diagonalization of the kinetic part (including the mass term) of the Lagrangian so as to permit usual perturbative analysis. While this would be a standard exercise, a slightly different formulation is more common Sakurai ; OConnell:1995nse , namely

[TABLE]

For the special point $g_{{}_{VJ}}=g_{{}_{V}}$ (necessary to maintain $F_{\rm em}(q^{2}=0)=1$), this transformation leads to, approximately,

[TABLE]

with the difference between the two Lagrangians being ${\cal O}(e^{3}/g_{{}_{V}}^{3})$. The existence, in eqn.(39), of a mass term for the photon is, of course, a concern. However, note that there now exists a mass-mixing between the photon and the $V$, in lieu of the earlier kinetic mixing. Summing to all orders in this mixing (no loops, though), the tree-level photon propagator can be easily seen to be proportional to

[TABLE]

While the $q^{2}\to 0$ limit reproduces the standard propagator (modulo a renormalization), the existence of the second pole is a reminder of the fact that the field redefinitions in eqn.(38) are neither complete nor even unitary. Neglecting this aspect for the time being, it is easy to see that, for $\gamma^{*}\to\pi^{+}\pi^{-}$, while the Lagrangian of eqn.(37) trivially reproduces $F_{\rm em}(q^{2}=0)=1$, the one in eqn.(39) would have done so only for $g_{{}_{VJ}}=g_{{}_{V}}$, as claimed earlier. With the absence of any direct coupling of the photon with the current $J^{\mu}$ existing in eqn.(39), the latter represents, explicitly, the situation of complete vector meson dominance, leading to a wide acceptance of such a phenomenological Lagrangian. Furthermore, in such a theory, $g_{{}_{VJ}}$ should be identical for all currents (involving fields of a given charge) so as to preserve gauge invariance.

With the second pole in the photon propagator appearing only at $m_{V}^{2}\,(1+e^{2}/g_{{}_{V}}^{2})$, the Lagrangians of eqns.(37&39) are expected to give identical results for $q^{2}\ll m_{V}^{2}$. Nonetheless, we will use the more common variant, namely eqn.(39). This can be trivially extended to include multiple vector mesons. Similarly, $J^{\mu}$ can be generalised to different currents. Given this, the amplitude for $\gamma^{*}\rightarrow\pi\rho$ can be expressed as

[TABLE]

such that

[TABLE]

In Table 2, we list the processes used to estimate the values of the different coupling constants. In calculating these, we adopt the natural leading order form for the three-point vector-vector-pseudoscalar meson vertex, viz. $\displaystyle\varepsilon_{\mu\nu\sigma\omega}p_{V1}^{\nu}p_{V2}^{\sigma}$.

Note that eqn.(41) gives a good approximation only for relatively small range of $q^{2}$-values where the two poles dominate in turn. Consequently, we use this form only for the region $q^{2}\leq 2.6\,{\rm GeV}^{2}$. For large $q^{2}$ values (q^{2}\mathrel{\hbox to0.0pt{\raise 2.15277pt\hbox{>}\hss}{\lower 2.15277pt\hbox{\sim}}}5\,{\rm GeV}^{2}), perturbative results (incorporating $k_{T}$ factorization etc.) are available and quite accurate (see Ref.Hua:2018kho ). In between these two regions, we interpolate, maintaining a $s^{-1}$ form (while poles do exist even in this region, their contributions are small on account of the corresponding vector-mesons having only suppressed couplings with a $\pi\rho$ pair). Finally, combining all the results, the electromagnetic form factor for $\pi\rho$, as used by us, is shown in Fig.6. It behoves us to estimate the accuracy of our calculation of the form factor. While comparison to data would be the natural check, unfortunately, there exists no experimental data for the $q^{2}>0.5\,{\rm GeV}^{2}$ region, thereby ruling out this possibility. For the $q^{2}<0.5\,{\rm GeV}^{2}$ region, data does exist, and a comparison with Ref.Maris:2002mz , does show some deviation, but never exceeding 25%. While even this may seem large, it should be realized that this particular range has very little effect on the numerical results that we present next. Rather, this should be considered as the maximal theoretical possible uncertainty in our calculations. Indeed, the errors are expected to be smaller for higher $q^{2}$ values, as the leading contributions have been accounted for.

5.4 Revisiting Relic abundance for vector interactions

We use these results to calculate the relic abundance and compare it with the case of free quarks(Fig.1). The contour satisfying $\Omega_{\varphi}h^{2}=0.1199$ is depicted in Fig.6. Once again, we show the physically untenable curve corresponding to the case of free quarks so as to facilitate easy comparison. Note that, for large $m_{\varphi}$, the realistic curve approaches the free-quark curve. This is quite similar to the case for the dimension-5 operators and is reflective of the fact that, for such $m_{\varphi}$ values, the free-quark approximation becomes better. More importantly, there exists a large peak with a substantial width just prior to the asymptotic region, where the free-quark approximation would have grossly underestimated the annihilation cross-sections and, hence, the sensitivity to the scale $\Lambda$. As can be realized from the behaviour of the form factor as displayed in Fig.6, this is but a reflection of the dominance of the $\phi$ and $\omega$ mesons. Despite the narrowness of the two resonances, their closeness implies that, when convoluted with the momentum spreads of the DM, they two peaks are no longer resolvable. Rather, the two contributions add coherently, with the first one dominating.

While the discussion above has concentrated only on the $\pi\rho$ final state, it is obvious that other final states too need to be taken into account for larger $m_{\varphi}$. Obvious candidates are states like $KK^{*}$, which dominate when the DM-pair couples to a strange-quark current. Relatable by $SU(3)$ symmetry (albeit broken badly) to the $\pi\rho$ final state, we can use the same formalism for this case too. More interesting are final states comprising two vector mesons, such as $\rho\rho$. A very similar analysis would go through for this as well. Indeed, the break in the fall of the dashed curve in Fig.6 and the subsequent rise for larger $m_{\varphi}$ owes its existence to the inclusion of the $\rho\rho$ state. The inclusion of even more states would drive this curve very close to the blue curve. This is only to be expected as, for $m_{\varphi}>2\,{\rm GeV}$, the annihilation can be well-approximated by quasi-free quarks.

6 CMB constraints

6.1 Effective relativistic degrees of freedom

Energy injection from DM annihilation in the early universe can alter the effective number of relativistic degrees of freedom $N_{\rm eff}$. Indeed, MeV-scale DM is especially constrained by these observations. If the DM freezes out after the neutrinos have decoupled (at $T=T_{\nu}^{\rm decoup}$), its annihilation will result in heating the $e^{-}$–$\gamma$ plasma relative to the neutrinos, thereby reducing the ratio of the neutrino and photon temperatures ($T_{\nu}/T_{\gamma}$). This results in a reduction of $N_{\rm eff}$ as $N_{\rm eff}\propto(T_{\nu}/T_{\gamma})^{4}$. From standard cosmology results, $N_{\rm eff}=3.046$ Ade:2015xua , and only small deviations from this value are allowed.

To find the expression for $T_{f}$ (equivalently, $x_{f}$), we can approximate $\langle\sigma\,v\rangle\sim\sigma_{0}(1+b/x_{f})$ (partial wave expansion of the cross section). On equating the interaction rate ($\Gamma(x_{f})$) with the expansion rate ($H(x_{f})$), we get

[TABLE]

Assuming that $Y\approx Y_{eq}^{\rm freeze-out}$ (or that $n_{\varphi}(T_{f})$ is equal to $n_{\varphi}$ today, i.e., at $T_{0}=2.73~{}{\rm K}$), we can solve this iteratively. For example, at the first iteration, the solution reads

[TABLE]

This implies that $x_{f}$ depends, mainly, on $m_{\varphi}$, $g_{\varphi}$ and $\sigma_{0}$. For a given particle, the magnitude of $\sigma_{0}$ is similar for all the operators under consideration, provided the respective Wilson coefficients are similar.

In Fig.7, we illustrate the decoupling temperature as a function of $m_{\varphi}$ for the $\mathcal{O}^{f}_{s}$ operator. With $\sigma_{0}$ being similar in magnitude for a given field, the value of $x_{f}$ would be very similar for the other operators too. As the figure shows, for $m_{\varphi}>20\,{\rm MeV}$, we have $T_{f}>$2.5 MeV. In other words, the scalar WIMP decouples prior to neutrino decoupling ($T_{\nu}^{\rm decoup}=2$ MeV). Consequently, on annihilation, it heats the neutrinos along with the photons and electrons, preserving the standard result of $N_{\rm eff}\approx 3.046$.

The situation seemingly gets complicated for $m_{\varphi}\in[6,20]$ MeV, when $T_{f}\approx T_{\nu}^{\rm decoup}$. However, note that the bulk of the entropy transfer due to DM annihilation still occurs at $T\sim m_{\varphi}/3$, i.e., prior to $\nu-$decoupling, and, hence, the model is safe from such constraints. Thus, for a complex scalar field $\varphi$, it is the $m_{\varphi}<$6 MeV range which is constrained by $N_{\rm eff}$. Similarly, for a real scalar, the limit is 3 MeV. These results are in consonance with those in Ref.Ho:2012ug . However, it should be reiterated that this is operative only in (standard) scenarios wherein the DM is presumed to have been produced thermally and having been in thermal equilibrium with the SM particles171717If light DM enters thermal equilibrium with the SM after neutrino-photon decoupling, then the constraints from measurements of $N_{\rm eff}$ are significantly relaxedBerlin:2017ftj ., for a long enough phase with its abundance being defined thereby. For a non-thermal DM, the couplings to neutrinos and/or photons will have to be tuned such that the model satisfies above constraints.

6.2 CMB observations and Indirect Detection

The Cosmic Microwave Background Radiation encodes information about the thermal history of the early universe, and is well described by SM physics. On the other hand, DM annihilation, at early times, into high energy photons or charged particles can not only heat the gas, but can also lead to atomic excitations and even its ionization. This increase in the amount of the ionized fraction causes an increase in the width of the last scattering surface, thereby affecting the power spectrum of the CMB Chen:2003gz ; Madhavacheril:2013cna .

The energy injected by DM annihilation into the CMBR depends on its number density $n_{\varphi}$ at that epoch, the rate of its annihilation into (charged) SM particles and the nature of the cascade of particles produced after annihilation. Due to this cascading, not the entire energy is transferred to the CMB (or the plasma in equilibrium with it) but only a fraction. To calculate the amount of the energy transferred, one needs to track the evolution of the hydrogen and helium ion fractions, and the spectra of $e^{\pm}$ and photons at that epoch. With these being temperature– and, hence, redshift–dependent, we are faced with a redshift dependent efficiency function $f(z)$ that describes the fraction of the energy absorbed by the CMB plasma. It has been argued Slatyer:2015jla , though, that the effect of $f(z)$ can be well-approximated by an effective, but redshift independent, efficiency function $f_{\rm eff}$. Indeed, Ref.2012PhRvD..85d3522F demonstrated that, given a set of $f(z)$ functions for a WIMP, the impact of an appropriately chosen $f_{\rm eff}$, on the CMB, is identical at the sub-percent level. This is the simplification that we shall adopt.

The rate of energy deposited, into the CMB, by DM pair annihilation per unit time per unit volume is given by:

[TABLE]

where $z$ is the redshift of the epoch, and $\rho_{c}$($\Omega_{\varphi}$) is the critical density of the universe (DM relic abundance) today, i.e., at $z=0$. The factor $(1+z)^{6}$ just encapsulates the standard evolution of the dark matter number density (note that the annihilation rate is proportional to $n_{\varphi}^{2}$). CMB observationsAde:2015xua constrain $P_{\rm ann}<4.1\times 10^{-28}$cm3s*-1GeV-1*. To translate this into the allowed region in the $m_{\varphi}-\Lambda$ plane, we need $f(z)$ (or, equivalently, $f_{\rm eff}$) and $\left<\sigma v\right>$. Here, $v$ is the relative velocity (in units of the velocity of light) of the second DM particle in the rest frame of the first. For a complex scalar field, the thermal average of the annihilation cross sections for different operators are given by

[TABLE]

For a real scalar field, similar expressions hold, but with an extra factor of $4$. As derived in appendix A, we have $\left<v^{2}\right>=6~{}T/m_{\varphi}$ and $\left<v^{4}\right>=60\,T^{2}/m_{\varphi}^{2}$.

The dependence of $\langle\sigma v\rangle$ on $\langle v^{n}\rangle$ is easy to understand in terms of the angular momenta, especially if the DM-pair is viewed as a composite pseudoparticle (with some angular momentum) decaying into a SM-pair. For $\mathcal{O}^{f}_{s}$ and $\mathcal{O}^{f}_{p}$, the amplitude has both s-wave and p-wave components, whereas for $\mathcal{O}^{f}_{v}$ and $\mathcal{O}^{f}_{a}$, no s-wave component can exist on account of the inherent angular momentum of the initial state. Similarly, for $\mathcal{O}_{\gamma}$ and $\mathcal{O}_{\tilde{\gamma}}$, no p-wave component may exist as it would require the diphoton state to exist in an antisymmetric state.

With the DM being nonrelativistic, the CMB constraints for p-wave annihilation are weaker compared to those for the cases driven by s-wave DM annihilation. For pure s-wave annihilation, $\left<\sigma v\right>$ is independent of velocity and hence $P_{\rm ann}$ is a redshift-independent parameter. For vector couplings, on the other hand, we need to estimate the velocity of DM at the epoch where these interactions are significant. To this end, we consider the epoch of kinetic decoupling. When the Hubble rate equates the rate of scattering, the dark matter can no longer reach kinetic equilibrium with the plasma through a high momentum exchange rate; the WIMPs kinetically decouple from the plasma and attain free-streaming. This allows us to write

[TABLE]

and use

[TABLE]

to calculate $\langle\sigma v\rangle_{\text{CMB}}$. As the DM was, hitherto, in kinetic equilibrium with the plasma, $\left<v^{2}\right>_{kd}=\sqrt{6\,T_{kd}/m_{\varphi}}$. A conservative estimate gives181818The exact value for different operators vary from $10^{-4}$ to $10^{-6}$, as can be estimated by equating the rate for elastic scattering (DM and SM) to the expansion rate, i.e., $n_{rel}\,\Gamma_{\rm elastic}\sim H(T_{d})$. $x_{kd}=T_{kd}/m_{\varphi}\sim 10^{-4}$. Using $z_{CMB}\sim 1100$, we may estimate $T^{\gamma}_{CMB}=T^{\gamma}_{\rm today}(1+z_{CMB})$. The temperature of $\varphi$ at recombination depends on the kinetic decoupling temperature. As long as the DM remains kinetically coupled to the plasma, we have $T^{\gamma}_{kd}=T^{\varphi}_{kd}$. Once the DM decouples kinetically, its temperature at a redshift $z\sim z_{CMB}$ is decided by its non-relativistic nature, i.e,

[TABLE]

This implies

[TABLE]

or, in other words, $\left<v^{2}\right>_{\text{CMB}}$ of $\varphi$ of $m=1\,{\rm MeV}$ exceeds that for $m=1\,{\rm GeV}$ by a factor of $10^{9}$.

As for the effective efficiency factor $f_{\rm eff}$, it depends upon the details of the model (in our case, the relative sizes of the Wilson coefficients), and rather than calculate it explicitly, we allow it to vary within the range $0.4<f_{\rm eff}<1$ which is commensurate with that advocated in Ref.Gonzalez-Morales:2017jkx . As we shall see later, our results are not going to be very sensitive to the exact choice.

Since DM annihilation to electrons and photons give us the tightest constraints, in Fig. 8, we depict the value of $\left<\sigma v\right>(\varphi\varphi\rightarrow e^{-}e^{+})$, as a function of $m_{\varphi}$ for the respective operators of eq.(1) and different final states. In Fig. 8, we depict the value of $\left<\sigma v\right>(\varphi\varphi\rightarrow e^{-}e^{+})$, as a function of $m_{\varphi}$ for the respective operators of eq.(1) and different final states. In each case, $\Lambda$ is chosen to be the maximum allowed for by the measurement of the relic density, viz. the condition $\Omega_{\varphi}\leq\Omega_{DM}$. Thus, it is the area above a curve that is allowed. Also depicted, in solid blue, is the curve corresponding to the aforementioned CMB observation viz. $P_{\rm ann}\sim 4.1\times 10^{-28}$cm3s*-1GeV-1*. The top (bottom) curves correspond to $f_{\rm eff}=0.4\,(1)$ respectively. Clearly, it is the area below these curves that is allowed. Had the quarks in the final state been truly free, we would, thus, have faced a seeming disagreement between the two sets of observations, at least for smaller values of $m_{\varphi}$. However, before we entirely discard such a mass range, we need to reconsider the correction wrought by considering bound states instead. As Fig.8 shows, this reduces the disagreement to a very large extent. The inclusion of even more bound states in the calculation of relic abundance, would have further reduced the remaining disagreement (especially for m_{\varphi}\mathrel{\hbox to0.0pt{\raise 2.15277pt\hbox{>}\hss}{\lower 2.15277pt\hbox{\sim}}}500\,{\rm MeV}). To appreciate this, recognize that such an inclusion would raise the DM annihilation cross section (into hadronic states) even further, thereby implying a raise in the preferred value of $\Lambda$ and, hence, a suppression in the $\varphi\varphi\to e^{+}e^{-}$ rates.

In addition, it should be borne in mind that the said disagreement depends not only upon our understanding of the early universe being perfect, but also on several key assumptions. For example, consider the case of the non-thermal DM, where the annihilation cross sections are very small and the final relic abundance is completely determined by its initial abundance which, in turn, depends on the model at hand. For such small cross sections, these limits can be evaded easily. For scenarios that fall somewhere in between the non-thermal and rigorously thermal DM, the constraints would need to be scaled appropriately.

7 Summary and Conclusions

In this work, we have systematically investigated the interactions of a light scalar DM particle with the SM sector within the framework of an effective field theory. This entails examining the constraints imposed by astrophysical/cosmological observations such as the CMB power spectrum, the deduced value of the relic abundance as well as the comparison of the photon and neutrino temperatures.

In the first part of the work, we considered leptons and free quarks as final states (while this is not a good approximation given the smallness of the DM mass and, hence, energies available, it serves to illustrate certain features and sets the stage for the second part of the study) and analysed the constraints imposed by the cosmologically deduced value of the relic density ($\Omega_{\rm DM}\,h^{2}$). This allowed us to obtain relatively robust upper limits on the scale $\Lambda$ of the effective field theory; as a higher value for $\Lambda$ would imply a smaller DM annihilation cross section, and, hence, a larger abundance. For the dimension-6 operators under consideration, we have $\Lambda\propto\sqrt{m_{\varphi}}$ so as to produce the right abundance, while, for the dimension-5 operators, $\Lambda\sim 10^{4}\,{\rm GeV}$ with only a weak dependence on $m_{\varphi}$. The value of $\Lambda_{\rm max}$ corresponding to the dimension-6 operators might, at the first sight, seem too low to have escaped detection in terrestrial experiments. However, this is not quite so as we would demonstrate in the accompanying paperpaper_2 .

As already stated, for DM masses below a couple of GeVs, it is not a good approximation to assume that a pair of DM particles may annihilate into a pair of free quarks as DM of such a mass would freeze out only around the QCD phase transition temperature, leaving them with very little overall energy. Consequently, we should, instead, re frame the analysis in terms of bound states. Given that baryon production is suppressed, we concentrate on mesons, devising appropriate methodology to determine the leading annihilation cross sections into such states. Chiral perturbation theory as well as techniques of dispersion analysis are used to obtain the effective couplings of DM with a host of light mesons, not only the pseudoscalars such as pions and kaons, but vectors as well. However, it is the pseudoscalar states that dominate for dimension-5 operators and obtaining the right relic abundance would need the scale of the effective theory to be related to the DM mass as $\Lambda\propto F(q^{2})m_{\varphi}^{-2}$. For dimension-6 operators, on the other hand, annihilation to final state (pseudo)scalar mesons is $v^{4}$ suppressed and, hence, we must include the vector mesons in the mix. For a final state comprised of a pair of vector mesons, very good results can be obtained using the analogues of the time-like electromagnetic form-factors. For a pseudo scalar-vector combination, on the other hand, a combination of data and the vector-meson-dominance model does the job. This leads to $\Lambda\propto\sqrt{m_{\varphi}\,F(q^{2})}$ for vector interactions. What is particularly heartening to see is that the inclusion of progressively more states brings the results closer and closer to that obtained with free quarks. This lends credence to the belief that the results found herein present a very good approximation and can be made even more robust by the inclusion of just a few more states at best.

An orthogonal constraint emanates from the requirement that the annihilation of the DM does not significantly alter the ratio of the neutrino and photon temperatures, an observable often recast as $N_{\rm eff}$ (or, the effective number of neutrino-like species). For the effective Lagrangians under consideration, once the requirement of reproducing the right relic abundance is imposed, the freeze-out temperature $T_{f}\sim d_{1}\,m_{\varphi}\,\exp(m_{\varphi}/d_{2})$ with the constants $d_{1,2}$ being only very weakly dependent on the exact nature of the current-current structure. A consequence is that the constraints are strong only for m_{\varphi}\mathrel{\hbox to0.0pt{\raise 2.15277pt\hbox{<}\hss}{\lower 2.15277pt\hbox{\sim}}}6\,{\rm MeV}, as for higher values of $m_{\varphi}$, the bulk of the entropy transfer to the plasma, takes place before $T\approx m_{\varphi}/3$, and, hence, before the neutrinos have decoupled. Even this constraint (for m_{\varphi}\mathrel{\hbox to0.0pt{\raise 2.15277pt\hbox{<}\hss}{\lower 2.15277pt\hbox{\sim}}}6\,{\rm MeV}) can be evaded if the DM had a non-thermal origin. However, with the dynamics of such DM being very sensitive to the spectrum and the structure of the theory, it does not easily lend itself to the effective Lagrangian treatment.

A competing constraint emanates from the shape of the CMB spectrum. The lack of significant distortions in the same puts an upper limit on the rate of DM annihilation ($P_{\rm ann}$) to, for example, $e^{\pm}$ or photon-pairs. The lower limit on $\Lambda$ that this translates to is, often, in ostensible contradiction with the aforementioned values of $\Lambda_{\rm max}$. These two opposing constraints, thus, seemingly rule out such a light DM (within the ambit of an effective theory). However, with the inclusion of bound states, the disagreement between the observables $\Omega_{\varphi}h^{2}$ and $P_{\rm ann}$ is rendered comparatively small. It should be realized, though, that existence of even such a small “discrepancy” depends crucially on the assumption that, in the early universe, the DM was in exact thermal equilibrium with the SM sector. If this restriction is released, or, in other words, a non-thermal initial condition on the DM allowed for, the constraints from $P_{\rm ann}$ are eased sufficiently enough to permit a large overlap with the parameter space allowed by $\Omega_{\rm DM}$.

Similar to the CMB constraints, the inclusion of bound states slightly changes the interpretation of the results in the direct detection experiments which we have discussed in the companion paper.

In summary, we may conclude that a substantial fraction of the parameter space of light scalar DM is still viable. Furthermore, an accurate estimation of the cosmological constraints needs the proper inclusion of bound states. While the inclusion of the few light mesons already given us a fast converging and robust result, the remaining uncertainties can be reduced in a straightforward (though painstaking) manner by the inclusion of even more states.

Appendix A Appendix A

For any observable $f(p_{1},p_{2},\dots,p_{n})$, constructed of the momenta of $n$ particles of a gas at equilibrium, the thermal average is given by

[TABLE]

where $E_{i}$ is the energy of the $i$’th particle. For a single non-relativistic species,

[TABLE]

Defining the relative velocity between two particles as

[TABLE]

we have

[TABLE]

Henceforth, we choose $k_{B}=1$. Effecting a change of variables, namely

[TABLE]

the Jacobian is unity, and

[TABLE]

Similarly,

[TABLE]