Dark matter from CP symmetry of order 4: evolution in the asymmetric regime

TL;DR

This paper investigates the thermal evolution of scalar dark matter candidates stabilized by a higher-order CP symmetry in a three-Higgs-doublet model, focusing on asymmetry, conversion processes, and potential indirect detection signals.

Contribution

It introduces a novel dark matter stabilization mechanism via CP symmetry of order four and analyzes its impact on dark matter evolution and observable signals.

Findings

Conversion processes affect dark matter asymmetry evolution.

Analytic and numerical solutions of Boltzmann equations are provided.

Potential for observable indirect detection signals is assessed.

Abstract

Multi-Higgs models equipped with global symmetries produce scalar dark matter (DM) candidates stabilized by the unbroken symmetry. It is remarkable that a conserved CP symmetry can also stabilize DM candidates, provided it is a CP symmetry of order higher than two. CP4 3HDM, the three-Higgs-doublet model with CP symmetry of order 4, is the simplest example of this kind. It contains two mass-degenerate scalar DM candidates $\varphi$ and $\bar\varphi$, each of them being a CP4 eigenstate and, therefore, its own antiparticle. A novel phenomenological feature of this model is the presence of $\varphi\varphi \leftrightarrow \bar\varphi\bar\varphi$ conversion process, which conserves CP. It offers a rare example of DM models in which self-interaction in the dark sector can significantly affect cosmological and astrophysical observables. Here, we explore the thermal evolution of these DM species in the asymmetric regime. We assume that a mechanism external to CP4 3HDM produces an initial imbalance of the densities of $\varphi$ and $\bar\varphi$. As the Universe cools down, we track the evolution of the asymmetry through different stages, and determine how the final asymmetry depends on the interplay between the conversion and annihilation $\varphi\bar\varphi \to $ SM and on the initial conditions. We begin with the analytic treatment of Boltzmann equations, present a detailed qualitative description of the process, and then corroborate it with numerical results obtained using a dedicated computer code. Finally, we check if the model can produce an observable indirect detection signal.

Figures

Equations

$$V$$

$$X=\left(\begin{array}[]{ccc}1&0&0\\ 0&0&i\\ 0&-i&0\end{array}\right)\,.$$

$$\phi_{1}=\left(\begin{array}[]{c}G^{+}\\ {1\over\sqrt{2}}(\mathrm{v}+h_{125}+iG^{0})\end{array}\right)\,,\quad\phi_{2}=\left(\begin{array}[]{c}H^{+}_{2}\\ {1\over\sqrt{2}}(H+ia)\end{array}\right)\,,\quad\phi_{3}=\left(\begin{array}[]{c}H^{+}_{3}\\ {1\over\sqrt{2}}(h+iA)\end{array}\right)\,,$$

$$m^2\equiv m_{h,a}^2$$

$$M^2\equiv m_{H,A}^2$$

$$m_{H^{\pm }}^2\equiv m_{H_2^{\pm },H_3^{\pm }}^2$$

$$\phi =\frac{1}{\sqrt{2}}(h+ia),\phi \stackrel{CP}{\to }i\phi .$$

$$\delta =\frac{n_{\phi }-n_{\bar{\phi }}}{n_{\phi }+n_{\bar{\phi }}},$$

$$\phi \stackrel{CP}{\to }i\phi ,{\phi }^{\ast }\stackrel{CP}{\to }-i{\phi }^{\ast },$$

$$\phi \phi ⟷\bar{\phi }\bar{\phi },$$

$$\{\begin{array}{l}\dot{Y}=-A(x)\left[Y\bar{Y}-Y_{\mathrm{eq}}^2(x)\right]-B(x)\left(Y^2-{\bar{Y}}^2\right)/4\\ \dot{\bar{Y}}=-A(x)\left[Y\bar{Y}-Y_{\mathrm{eq}}^2(x)\right]+B(x)\left(Y^2-{\bar{Y}}^2\right)/4,\end{array}$$

$$A(x)=\sqrt{\frac{g(x)\pi }{45}}\frac{m{M}_{\mathrm{Pl}}}{x^2}⟨\sigma v{⟩}_{\mathrm{SM}}(x),$$

$$B(x)=\sqrt{\frac{32g(x)}{45}}\frac{m{M}_{\mathrm{Pl}}}{x^{5/2}}{\sigma }_{\mathrm{conv}}.$$

$$A(x)\approx \frac{a}{x^2},B(x)\approx \frac{b}{x^{5/2}},$$

$$Y_t=\frac{Y+\bar{Y}}{2},\Delta =\frac{Y-\bar{Y}}{2},\delta =\frac{\Delta }{Y_t}.$$

$$\{\begin{array}{l}\dot{Y_t}=-A(x)\left[Y_t^2-{\Delta }^2-Y_{\mathrm{eq}}^2(x)\right]\\ \dot{\Delta }=-B(x)Y_t\Delta .\end{array}$$

$$\{\begin{array}{l}\dot{Y_t}=-A(x)\left[Y_t^2(1-{\delta }^2)-Y_{\mathrm{eq}}^2(x)\right]\\ \dot{\delta }=-\delta \left[B(x)Y_t+\dot{Y_t}/Y_t\right].\end{array}$$

$$\dot{\delta }=\delta \cdot Y_t\left[A(x)\left(1-{\delta }^2-\frac{Y_{\mathrm{eq}}^2}{Y_t^2}\right)-B(x)\right].$$

$$Y_{\mathrm{eq}}(x)=c(x)x^{3/2}\exp (-x),c(x)=\frac{90}{g(x)(2\pi {)}^{7/2}}.$$

$$\dot{\delta }\approx \delta \left(1-\frac{3}{2x}\right).$$

$$\delta (x)\propto \exp (x).$$

$$Y_t(x\to \infty )=Y_t(x_f){\left(1+a\frac{Y_t(x_f)}{x_f}\right)}^{-1}\approx \frac{x_f}{a}.$$

$$\dot{\delta }\approx \delta \left(1-{\delta }^2\right).$$

$$\delta (x)\approx \frac{{\delta }_0\exp (x-x_0)}{\sqrt{1+{\delta }_0^2\exp (2x-2x_0)}},$$

$$x_f\approx 1+\ln \left(\frac{7}{{\delta }_0}\right),$$

$$r\equiv \frac{\bar{Y}}{Y+\bar{Y}}=\frac{1-\delta }{2}.$$

$$\dot{r}=-2\Delta A(x)\left(r-\frac{Y_{\mathrm{eq}}^2}{4{\Delta }^2}\right).$$

$$r(\infty )\approx \frac{\sqrt{\pi }}{2}{a}^2{c}^2(a\Delta {)}^{-1/4}\exp (-4\sqrt{a\Delta }),$$

$$\delta(x)\approx\delta_{0}\exp\left[(x-x_{0})-m\tilde{M}_{\rm Pl}\sigma_{\rm conv}\Big{(}{\rm E_{1}}(x_{0})-{\rm E_{1}}(x)\Big{)}-\frac{{\displaystyle 3}}{{\displaystyle 2}}\ln(x/x_{0})\right]\,.$$

$$E_1(x)={\int }_x^{\infty }\frac{dt}{t}{e}^{-t},$$

$$m{\tilde{M}}_{\mathrm{Pl}}{\sigma }_{\mathrm{conv}}{e}^{-x_{\ast }}=x_{\ast }.$$

Full text

I Introduction

I.1 Dark matter candidates from the Higgs sector

Many astrophysical and cosmological observations offer a compelling evidence of the existence of Dark Matter (DM) Hinshaw:2012aka ; Ade:2015xua . From the particle physics perspective on DM, the main paradigm is to assume that the DM candidate is stable on cosmological time scales, cold, i.e. non-relativistic at the onset of galaxy formation, non-baryonic, neutral, and weakly interacting, and that its relic abundance is obtained via thermal freeze-out in the early Universe, see e.g. Bertone:2004pz for a review. Such particles are absent in the Standard Model (SM). Together with the absence of an equally compelling evidence for scattering of DM on nuclei Akerib:2016vxi ; Aprile:2017iyp or for astrophysical signals which would unequivocally confirm DM annihilation processes TheFermi-LAT:2017vmf , it motivates the vast field of building models beyond the SM which include DM candidates.

One popular class of such models is represented by symmetry-based multi-Higgs models, in which DM candidates arise naturally Ivanov:2017dad ; Ilnicka:2018def . If a model features inert scalars — that is, fields not involved in generation of masses of gauge bosons or fermions — and if the inert scalar sector transforms non-trivially under a new global symmetry, which remains unbroken upon minimization of the scalar potential, then the lightest inert scalar is automatically protected against decay and serves as a DM candidate. One of the simplest and actively studied examples is the Inert Doublet Model (IDM), the version of the two-Higgs Doublet Model (2HDM) with the exact $\mathbb{Z}_{2}$ global symmetry Deshpande:1977rw ; Cao:2007rm ; Barbieri:2006dq ; LopezHonorez:2006gr which leads to interesting astrophysical and collider signatures Ilnicka:2015jba ; Belyaev:2016lok ; Belyaev:2018ext ; Kalinowski:2018ylg ; Kalinowski:2018kdn .

Scalar DM models with more than two Higgs doublets have also attracted much attention. In particular, within three-Higgs-doublets models (3HDMs), one can incorporate several global symmetry groups Ivanov:2012ry ; Ivanov:2012fp , which lead to a variety of options in DM model building Keus:2014jha ; Ahriche:2015mea ; Keus:2015xya ; Cordero-Cid:2016krd . It is interesting to remark that in almost all multi-Higgs models it is extremely difficult to probe and constrain self-interaction within the dark sector. It can be rather strong and can bear certain patterns, but usually it does not affect the DM relic density nor can it be probed through the direct or indirect detection experiments, see for example a detailed analysis in the case of IDM Sokolowska:2011sb ; Sokolowska:2011aa .

I.2 Dark matter stabilized by an exotic $CP$ symmetry

In Ivanov:2015mwl , another 3HDM was described, in which the only symmetry imposed was the $CP$ symmetry of order 4 (CP4). It is the first model featuring such a novel, physically distinct Haber:2018iwr form of $CP$ symmetry without producing any accidental symmetries. If CP4 remains unbroken after minimization of the Higgs potential, the model leads to a pair of mass-degenerate scalar DM candidates, which are stabilized, remarkably, by the $CP$ invariance in the scalar sector. This non-standard stabilization mechanism may have peculiar consequences for the DM evolution, and this is what we explore in the present work.

Let us briefly describe the setting of this model, leaving subtle details for the Appendix. CP4 3HDM is a three-Higgs-doublet model based on Higgs doublets $\phi_{1,2,3}$ equipped with a generalized $CP$ symmetry of order 4 denoted CP4 Ivanov:2011ae ; Ivanov:2015mwl . There are no new fields apart from the new Higgs doublets. Its potential can be written as

[TABLE]

Here, all parameters apart from $\lambda_{8,9}$ are real. This potential is invariant under the CP4 transformation $\phi_{i}(\vec{r},t)\mapsto X_{ij}\phi_{j}^{*}(-\vec{r},t)$, where matrix $X$ is

[TABLE]

The potential $V$ is the most general potential compatible with CP4, and it represents the minimal model which realizes this novel form of $CP$ symmetry.

Although multi-Higgs-doublet models offer a large basis change freedom, there exists no basis choice in which CP4 would take the “standard” form $\phi_{i}\mapsto\phi_{i}^{*}$. Or equivalently, there is no basis in which all parameters of the potential would become real. Technically, this is so because applying CP4 twice results not in the identity transformation but in a family symmetry $XX^{*}=\mathrm{diag}(1,\,-1,\,-1)$. One needs to apply CP4 four times to get the identity transformation. Thus, CP4 is a genuinely novel form of $CP$ symmetry, physically distinct from the “standard” $CP$. This distinction was also established at the level of observables Haber:2018iwr .

In a detailed scan of the parameter space of this model, one would need to satisfy the boundedness-from-below conditions and, to make the perturbative calculation reliable, the perturbative unitarity constraints. The necessary and sufficient boundedness-from-below conditions are not yet known, but a simple set of sufficient conditions was proposed in Ref. Ferreira:2017tvy . The perturbative unitarity constraints for CP4 3HDM were written and studied in Bento:2017eti ; Kopke:2018vyw . However, since we do not perform a detailed scan in this work, we will not use these conditions.

We are interested in the regime where CP4 is preserved in the vacuum and protects the scalar DM candidates against decay. The only vacuum expectation value (vev) alignment which conserves CP4 is $\langle\phi_{1}^{0},\phi_{2}^{0},\phi_{3}^{0}\rangle=(\mathrm{v},0,0)/\sqrt{2}$. The vev symbol $\mathrm{v}=246$ GeV is not to be confused with the velocity $v$. We expand the doublets around this vacuum as

[TABLE]

where $G^{+}$ and $G^{0}$ represent the would-be Goldstone bosons which are absorbed into the gauge bosons $W^{+}$ and $Z$ by the Brought-Englert-Higgs mechanism.

The SM-like Higgs boson comes from the first doublet and has mass $m_{h_{125}}^{2}=2\lambda_{1}\mathrm{v}^{2}=2m_{11}^{2}$. All other scalars coming from the two inert doublets show a pairwise degenerate spectrum:

[TABLE]

where we introduced the shorthand notation $\lambda_{346}=\lambda_{3}+\lambda_{4}-\lambda_{6}$ and $\bar{\lambda}_{346}=\lambda_{3}+\lambda_{4}+\lambda_{6}$. Notice that the pair-degenerate mass spectrum is imposed by a discrete, not continuous symmetry Haber:2018iwr . Due to this degeneracy, the spectrum of the model resembles a typical 2HDM spectrum: three values of neutral Higgs masses and one value of charged Higgs mass. However, the counting of additional scalars is different from 2HDM, which can be detected in principle by different production and decay rates, see more discussion in Haber:2018iwr .

The two neutral fields $h$ and $a$ are the DM candidates protected by CP4 symmetry. The relation $m<M$ implies $\lambda_{6}>0$ so that $\lambda_{346}<\bar{\lambda}_{346}$. Also, we assume that $\lambda_{6}-\lambda_{4}>0$ to guarantee that the lightest inert scalars are neutral. The heavier inert scalars become unstable decaying into the lightest ones and into the SM fields. They do not directly affect the relic density evolution.

Although $CP$ is conserved, one cannot classify the neutral scalars $h$ and $a$ as $CP$-even or $CP$-odd. They are not even $CP$-eigenstates, as they transform under CP4 as $h\xrightarrow{CP}-a$ and $a\xrightarrow{CP}h$. However, one can combine them into CP4 eigenstates Ivanov:2015mwl ; Aranda:2016qmp :

[TABLE]

The fields $\varphi$ and $\varphi^{*}$ are $CP$-eigenstates but they are characterized by a $CP$-charge defined modulo 4 rather than $CP$-parity. In order for this construction to be valid, the fields $h$ and $a$ cannot possess any other conserved quantum numbers, see Aranda:2016qmp for a detailed discussion.

I.3 CP4 3HDM in the asymmetric DM regime

A first discussion of the qualitative features of the DM sector in CP4 3HDM was given in Ivanov:2015mwl ; Haber:2018iwr . A more detailed quantitative insight into the DM evolution was given in Kopke:2018vyw . In this paper, we present a detailed exploration of another aspect of DM in CP4 3HDM — its evolution in the asymmetric DM regime.

The thermal evolution of asymmetric DM has been the subject of numerous studies, see Graesser:2011wi ; Iminniyaz:2011yp ; Zurek:2013wia ; Petraki:2013wwa and references therein. The main idea is that, at some point in the thermal history of the Universe, an asymmetry between DM particles and antiparticles was created, probably along with the baryon asymmetry of the Universe. Then one has two main options. If the DM particles $\varphi$ are associated with a new conserved global or gauge $U(1)$ charge, then the asymmetry $n_{\varphi}-n_{\bar{\varphi}}$ remains unchanged in the later evolution, which typically leads to an almost complete domination of DM particles over antiparticles. In terms of the relative asymmetry $\delta$ defined as

[TABLE]

this regime corresponds to $\delta\to 1$ at present. A typical feature of these models is a very suppressed or absent present-day annihilation signal. If the DM antiparticles have already died out, the DM particles cannot annihilate anymore. Another possibility is that there exists a novel interaction which violates the DM number conservation. If such violation is present in the quadratic terms, it will induce oscillations between DM particles and antiparticles and will tend to wash out surviving asymmetry, unless the oscillations are extremely weak Buckley:2011ye ; Tulin:2012re ; Cirelli:2011ac .

We show in this paper that the asymmetric DM evolution in CP4 3HDM has characteristic features which lead to a richer list of options. The lightest inert particles in CP4 3HDM are $\varphi\equiv\varphi^{*}|0\rangle$ and $\bar{\varphi}\equiv\varphi|0\rangle$, which have equal masses but which are not particle and antiparticle of each other. Each of them is an eigenstate of the $CP$ symmetry, since the conserved symmetry CP4 acts on the complex fields $\varphi$ and $\varphi^{*}$ as

[TABLE]

see more details in the Appendix. Thus, the asymmetric DM regime in CP4 3HDM is not directly related to the particle-antiparticle asymmetry. It is rather an asymmetry between two sorts of mass-degenerate complex fields.

Having said that, we admit that it is easier to discuss the properties of $\varphi$ and $\bar{\varphi}$ by calling them particle and antiparticles, respectively. This is a language abuse, but we will stick to it in this paper. We stress once again that, despite using this language, we recognize that $\varphi$ and $\bar{\varphi}$ are their own antiparticles.

The transformation law under CP4 given by Eq. (7) allows us to label the two particles not with $CP$ parities but with $CP$ charges $q$ defined modulo 4. This $CP$ charge protects them against decay and against $\varphi\leftrightarrow\bar{\varphi}$, as well as against annihilation $\varphi\varphi\to Z$.111If the heavier inert states are taken into account, there appears a number of coannihilation channels. However, if these scalars are significantly heavier, they quickly decay and do not affect the thermal freeze-out process nor the late-time evolution. However, the scalar potential contains terms that induce the conversion process (see the Appendix for the conversion cross section):

[TABLE]

which is well compatible with $CP$ charge conservation: $+2$ is equal to $-2$ modulo 4. Therefore, although the DM oscillations are forbidden, $\bar{\varphi}$’s can still be regenerated in $\varphi$-dominated Universe, pick up a $\varphi$ and annihilate into SM fields. Notice that this conversion is governed exclusively by the quartic self-interaction in the DM sector.

The relative intensity of the conversion and annihilation processes will affect the present day relic asymmetry as well as the total DM abundance. Depending on the parameters of the model and on the initial asymmetry222We do not discuss the origin of this primordial asymmetry and its possible relation with the baryon asymmetry of the Universe. It certainly goes beyond CP4 3HDM per se and requires additional phenomena, which may exist in a more comprehensive model at higher energy and leads to CP4 3HDM at the electroweak scale. In this work, we simply postulate an initial asymmetry at certain temperature and then investigate its evolution during thermal freeze-out., the outcome of the thermal evolution may be almost fully symmetric, almost total asymmetric, or partially asymmetric. We will trace the boundary between the extreme regimes, link the final asymmetry with its primordial value, and investigate the prospects of the indirect detection signal in these regimes. Thus, this model offers a rare example where self-interaction in the DM sector does significantly affect the cosmological observables.

The structure of the paper is as follows. In Section II we write down Boltzmann equations and discuss the qualitative features of their solutions. We show that, starting from the relatively low values of $x\equiv m/T\sim 1$, the relative asymmetry $\delta$ first drops (depletion stage), then steadily grows (recovery stage), and either approaches unity or flattens at some value. In Section III, we corroborate these results with the numerical analysis performed by means of a dedicated evolution code. Then, in Section IV, we estimate the annihilation signal in various regimes, and finally we draw our conclusions. The Appendix provides some details of CP4 3HDM and the properties of its DM candidates.

II Asymmetry evolution: qualitative analysis

II.1 Boltzmann equations

The evolution of particles $\varphi$ and $\bar{\varphi}$ (which, in abuse of language, will often be called below particle and antiparticle) in the homogeneous and isotropic expanding Universe is described by a system of coupled Boltzmann equations

[TABLE]

where the factor $1/4$ is introduced for future convenience. In these equations, $Y$ and $\bar{Y}$ denote the number density in the comoving volume of $\varphi$ and $\bar{\varphi}$, respectively. For the sake of brevity, from now on we refer to these quantities as “densities”, except for Section IV, where we will need to deal with physical number densities. The dot above a symbol denotes the derivative with respect to $x$, e.g. $\dot{Y}\equiv dY/dx$, where $x=m/T$, $m$ is the mass of $\varphi$ and $T$ is the temperature of the plasma. Eqs. (9) explicitly mention the $x$ dependence of the parameters; of course, $Y$ and $\bar{Y}$, being variables, also depend on $x$, which is implicitly assumed in all equations. $Y_{\rm eq}(x)$ is the equilibrium density at given $x$ of both $\varphi$ and $\bar{\varphi}$ in the $\varphi/\bar{\varphi}$ symmetric plasma. In an asymmetric situation, the actual equilibrium densities of $\varphi$ and $\bar{\varphi}$ differ but their product is still equal to $Y_{\rm eq}^{2}$. Finally, $A(x)$ and $B(x)/4$ denote the intensities of the annihilation $\varphi\bar{\varphi}\leftrightarrow\,{\rm SM}$ and conversion $\varphi\varphi\leftrightarrow\bar{\varphi}\bar{\varphi}$ processes, respectively. The first one is determined as follows333The details on the derivation of the exact velocity-averaged cross sections and equilibrium number density can be found, for instance, in the Appendix of Ahmed:2017dbb .

[TABLE]

where $g$ denotes the number of relativistic degrees of freedom at the moment $x$, $M_{\rm Pl}$ is the Planck mass, and $\langle\sigma v\rangle_{\!\,{\rm SM}}$ is the total velocity-averaged cross section of all interactions that turn $\varphi$ and $\bar{\varphi}$ into SM particles. A similar expression holds for $B(x)$, where one uses the corresponding velocity-averaged cross section $\langle\sigma v\rangle_{\rm conv}$ instead. Notice that $A(x)$ and $B(x)$ are expected to show distinct $x$ dependences. In the case of $s$-wave annihilation of slow DM particles into lighter SM fields, $\langle\sigma v\rangle_{\!\,{\rm SM}}(x)\approx$ const. For the $\varphi\varphi\leftrightarrow\bar{\varphi}\bar{\varphi}$ conversion process, due to $m_{\varphi}=m_{\bar{\varphi}}$, it is the cross section $\sigma_{\rm conv}$ that is approximately constant, see Eq. (51). Thus, taking into account the relation between the average velocity and temperature Gondolo:1990dk ; Cannoni:2013bza , the velocity-averaged cross section is $\langle\sigma v\rangle_{\rm conv}=\sigma_{\rm conv}\langle v\rangle=\sigma_{\rm conv}\cdot 2/\sqrt{\pi x}$. Finally, we arrive at the following expression

[TABLE]

Neglecting the slow change of $g$ with $x$, one can get the following approximation

[TABLE]

which shows that the relative importance of the conversion process to annihilation to SM fields decreases as $x$ grows.

Let us now split the two densities into the symmetric and asymmetric contributions:

[TABLE]

Here, $Y_{t}$ stands for half the total density of $\varphi$ and $\bar{\varphi}$, which we will call for brevity just the total density, $\Delta$ stands for half their difference, while $\delta$ quantifies the relative amount of asymmetry between $\varphi$ and $\bar{\varphi}$. The fully symmetric case corresponds to $\delta=0$, while the complete asymmetry implies $\delta=1$ or $-1$. For concreteness, we will assume that $\varphi$ dominate over $\bar{\varphi}$ so that the asymmetry is always non-negative.

By taking the sum and the difference in (9), one can recast these equations into evolution of $Y_{t}$ and $\Delta$:

[TABLE]

One sees that the conversion process washes out the absolute value of the asymmetry. However, since the total density also drops, the evolution of the relative asymmetry may be quite sophisticated. To describe it, we further recast (14) into equations for $Y_{t}$ and $\delta$:

[TABLE]

The last equation can also be explicitly expanded as

[TABLE]

One observes that there are several competing terms inside the brackets of (16). Thus, the relative asymmetry can grow, decrease, or change the regime during the $x$ evolution. It is this evolution that we will investigate below.

II.2 Qualitative analysis without conversion

Before presenting numerical results, let us discuss the expected qualitative features of the $x$ evolution of the relative asymmetry. We begin with the case when no conversion takes place, $B(x)=0$. In this case the evolution of the relative asymmetry $\delta$ can still be non-trivial, since the bracket in Eq. (16) can, in principle, change sign. Although the main features of the evolution are well known in this case, see e.g. Graesser:2011wi , we outline them to set up the notation.

Suppose one starts at $x=x_{0}={\cal O}(1)$ with a non-zero initial $\delta_{0}$ and some initial density $Y_{t0}$. These initial conditions, $\delta_{0}$ and $Y_{t0}$, must be considered as free parameters, since the origin and the amount of the primordial asymmetry go beyond CP4 3HDM and need to be derived from a theory at a higher mass scale, which simplifies to CP4 3HDM at the electroweak scale.

At first, when $x$ is well below the freeze-out point for $Y_{t}$, the total density $Y_{t}$ quickly approaches the equilibrium density $Y_{\rm eq}/\sqrt{1-\delta^{2}}$, either from above or from below, depending on $\delta_{0}$ and $Y_{t0}$. In particular, if we start with $\delta_{0}$ close to 1 and $Y_{t0}\sim Y_{\rm eq}$, we can see a brief drop in $\delta$ due to a momentary increase of $Y_{t}$. This equilibration episode, which is driven by the choice of the initial conditions, will be quickly taken over by the steady growth of $\delta$.

The subsequent evolution of $\delta$ can be deduced directly from (15). If $\delta_{0}\ll 1$ and $\delta$ stays small during thermal evolution, then, well before the freeze out, $Y_{t}\approx Y_{\rm eq}$ and is approximately determined by the Maxwell-Boltzmann distribution:

[TABLE]

Plugging this into the second equation in (15) one immediately gets

[TABLE]

For values of $x\gtrsim 10$, it further simplifies to $\dot{\delta}\approx\delta$ and results in a very simple dependence:

[TABLE]

If $\delta$ still stays small during this entire evolution until the freeze-out point for the total density $x=x_{f}$, one obtains the growth factor of the order of $\exp(x_{f}-x_{0})$.

After the freeze-out for the total density, $Y_{t}$ reaches a plateau, and $Y_{\rm eq}/Y_{t}$ quickly drops and can be neglected. If $\delta$ is still considerably less than $1$, the first equation in (15) can be explicitly integrated and one gets

[TABLE]

In the last passage here, we took into account that the extra term in the brackets is typically ${\cal O}(10)\gg 1$. Just like in the usual symmetric case, the freeze-out moment $x_{f}\sim 25$ is set by the observed relic density. Since $\Delta=\mathrm{const}$ in the absence of conversion, we have $\delta\propto Y_{t}^{-1}$, and the late-time value of the asymmetry $\delta_{\infty}$ grows from $\delta(x_{f})$ by the same factor as in Eq. (20). Thus, it is not guaranteed that $\delta$ eventually reaches unity: if the initial value $\delta_{0}$ was sufficiently small, $\delta$ can stay small throughout the entire freeze-out process.

It may also happen that $\delta$, being initially moderately small, grows close to unity during the above evolution. In this case, the $1-\delta^{2}$ factor in (15) changes the dynamics of $Y_{t}$ and shifts the freeze-out moment $x_{f}$ to earlier times. One can estimate the moment when $\delta$ approaches unity and the total density effectively freezes out. For that, we restore the non-zero $\delta$ and, during the thermal equilibrium stage, write the total density as $Y_{\rm eq}/\sqrt{1-\delta^{2}}$. The ratio $\dot{Y_{t}}/Y_{t}$ will now contain contributions not only from $Y_{\rm eq}$ but also from $\delta$. Substituting it again in the second equation of (15) and omitting the $3/2x$ term for simplicity, we obtain after some algebra:

[TABLE]

Integrating the above equation, we get

[TABLE]

which shows how $\delta(x)$ approaches unity.

For a rough estimation, let us define the freeze-out moment for the total density as the moment at which $\delta=0.99$. We then obtain

[TABLE]

where for definiteness we set the initial moment $x_{0}=1$. We see that any sizable initial asymmetry $\delta_{0}>10^{-8}$ shifts the total density freeze out to $x_{f}\lesssim 22$, which is earlier than in the symmetric case.

Notice that the moment $x=x_{f}$ only settles $Y_{t}\approx\Delta$ in the case of the early freeze-out, while the density of antiparticles and, therefore, $\delta$ keep evolving. In this regime, the relevant variable is not $\delta$ itself, but the relative amount of antiparticles:

[TABLE]

For $x>x_{f}$, it still keeps decreasing exponentially, $r\approx 5\cdot 10^{-3}\exp[-2(x-x_{f})]$, and it settles to $r(\infty)$ only when the evolution departs from thermal equilibrium. The evolution equation for $r\ll 1$ takes the following form

[TABLE]

Solving it and analyzing the solution with Laplace’s method keeping in mind that $a\Delta\gg 1$, one finds that $r$ begins flattening out around $x\sim x_{r}=\sqrt{a\Delta}$, and its asymptotic value is exponentially suppressed as

[TABLE]

where $c$ is the prefactor in the equilibrium density, Eq. (17).

It is important to appreciate that the value of $Y_{t}(\infty)\approx\Delta$ and the parameter $a$ are not related anymore as it was in Eq. (20). One can, for example, consider the annihilation cross section parameter $a$ much larger than what is typically needed for a given mass in the symmetric DM scenario for the same mass. The final total density will still be very close to $2\Delta$, but the relative contribution of antiparticles will be strongly suppressed according to (26).

II.3 Switching on conversion

Let us now include the conversion process, $B(x)\not=0$. Its role is to wash out asymmetry, and the question is whether this washing out is sufficiently strong to modify the dynamics of the asymmetry. Since the conversion process is determined by a combination of quartic self-coupling constants in the dark sector (see Appendix) this situation is the rare example of DM model where interactions in the dark sector have strong impact on the cosmological observables.

The first observation is that the conversion process may lead to a pronounced depletion of the asymmetry at the early stage of evolution. The entire thermal history of the asymmetry may now have three distinct stages, schematically shown in Fig. 1. At the beginning of the evolution, the initial value of $\delta$ drops down due to the conversion (depletion stage). Then, while the density and the conversion cross section decrease, the conversion process becomes ineffective and $\delta$ reaches its minimum at $x_{*}$ and starts to grow (recovery stage). This is because the equilibrium symmetric part of the total density decreases due to the expansion. After the total density freezes out the proportion of particles and antiparticles ceases to change as well and $\delta$ approaches its constant value (flattening stage). Note, that the beginning of the recovery stage also corresponds to the moment, when the value of the absolute asymmetry $\Delta$ becomes constant.

If, during the depletion stage, the two processes are of the same intensity, $B\sim A$, the $B(x)$ term in (16) strongly dominates and forces $\delta$ to sharply decrease all the way to the freeze-out moment, skipping the recovery stage and leading to a vanishingly small final asymmetry. If one wants the asymmetry to recover, one must assume that $B$ is not only smaller, but significantly smaller than $A$. Indeed, during the depletion stage, the coefficient $A$ stands in front of $1-\delta^{2}-Y_{\rm eq}^{2}/Y_{t}^{2}\ll 1$, and the conversion process competes with this strongly reduced annihilation. As $x$ grows, $\delta$ decreases and $Y_{t}$ starts to depart from $Y_{\rm eq}$. The coefficient $1-\delta^{2}-Y_{\rm eq}^{2}/Y_{t}^{2}$ steadily grows until the moment $x=x_{*}$ when the bracket in (16) changes sign, $\delta$ reaches its minimum value $\delta_{*}$, and the recovery stage begins.

As long as $\delta$ stays small $Y_{t}\approx Y_{\rm eq}$, one can explicitly integrate the second equation in (15) to obtain the following behavior:

[TABLE]

To keep the expressions compact, we introduce here $\tilde{M}_{\rm Pl}=\sqrt{45/g\pi^{7}}\,M_{\rm Pl}$, which differs from the commonly used $M^{*}_{\rm Pl}$ defined through the relation $H(T)=M^{*}_{\rm Pl}/T^{2}$ as $\tilde{M}_{\rm Pl}=2M^{*}_{\rm Pl}/\pi^{2}$. Assuming that $g\sim 100$ and does not change substantially during the thermal evolution, we estimate $\tilde{M}_{\rm Pl}\approx 10^{-2}M_{\rm Pl}\approx 10^{17}$ GeV. Finally, ${\rm E_{1}}(x)$ is the exponential integral defined as

[TABLE]

which can be approximated by $E_{1}(x)\approx e^{-x}/x$ already for moderately large values of $x$.

The depletion and recovery stages arise from the interplay of the first two terms in the exponential of Eq. (27). At moderate values of $x$, the sharply decreasing second term drives $\delta$ down. We highlight here the fact that the magnitude of the drop in $\delta$ is exponentially sensitive to the value of the conversion cross section. When this drop flattens, the growing first term takes over and leads, in the recovery stage, to a steady growth of $\delta$ similar to Eq. (19).

The position of the minimum $x_{*}$ can be read off Eq. (27). Neglecting the slowly varying $\ln(x/x_{0})$ term, one obtains a transcendental equation for $x_{*}$:

[TABLE]

To the logarithmic accuracy, we can estimate the solution as

[TABLE]

We could get the same expression just by looking at the moment when the absolute value of the asymmetry $\Delta$ freezes out, analogously to the common freeze out of the DM density in the symmetric case444See Section 9.2 in gorbunov2011introduction for the derivation of the formula, similar to Eq. (30).. Indeed, during the recovery stage, the conversion process effectively switches off, and it happens slightly after $x_{*}$. The value of $\Delta$ at this point can be estimated as

[TABLE]

where we neglected ${\rm E_{1}}(x_{*})\ll{\rm E_{1}}(x_{0})$ for $x_{*}\gtrsim 3$.

The subsequent evolution depends on the value of $\delta$ at this moment. If $\delta$ stays small, one expects approximately the same exponential grows of $\delta$ as in Eq. (19). At $x=x_{f}$, the moment of freeze-out for the total density, $\delta$ increases with respect to $\delta_{*}$ by the factor $\exp(x_{f}-x_{*})$. If $\delta$ still remains small at this moment, the flattening stage unfolds, during which $\delta$ experiences an additional residual increase by the same factor as in Eq. (20). Alternatively, if $\delta$ approaches unity, then the antiparticle contribution dies out similarly to the previous case, with $\Delta$ being effectively frozen at $x\approx x_{*}$.

II.4 Locating the boundary

We showed above that there are two extreme regimes: weakly asymmetric, with the final $\delta\ll 1$, and strongly asymmetric, when $\delta\approx 1$ and $r\ll 1$ at the end of thermal evolution. Let us now locate the boundary between the two regimes.

We see that the stronger is the conversion cross section $\sigma_{\rm conv}$, the more pronounced is the drop and the later asymmetry begins to restore. If this happens rather close to the total density freeze out, the asymmetry does not have enough time to recover to any considerable value. One can estimate the value of $\sigma_{\rm conv}$ which can lead, starting from a small initial $\delta_{0}$, to a considerable asymmetry around the total density freeze out. One just needs to extend Eq. (27) all the way to $x_{f}$ and require that the final result be of order 1:

[TABLE]

Notice that we neglected here $E_{1}(x_{f})$ with respect to $E_{1}(x_{0})$. For $x_{0}=1$, $\delta_{0}=0.1$ and $x_{f}\sim 20\div 30$ we numerically obtain $x_{*}\sim 3$. It confirms that the conversion interaction should freeze out relatively early compared to the total density freeze out, in order for a sufficient asymmetry to survive at the end.

Next, using Eq. (29), we can recast this condition into an upper bound of the conversion cross section:

[TABLE]

For the value of $x_{*}$, obtained above, and for $m\sim 100\,{\rm GeV}$, this relation gives $\sigma_{\rm conv}\lesssim 10^{-17}\,{\rm GeV}^{-2}$. Comparing this result to the typical velocity-averaged cross section of the interaction between WIMPs and SM particles needed to reproduce the correct relic abundance $\sigma_{\rm SM,0}=2.5\cdot 10^{-9}\,{\rm GeV}^{-2}$, one sees that the conversion cross section has to be about eight orders of magnitude smaller. Within CP4 3HDM, the cross section of the conversion process $\varphi\varphi\leftrightarrow\bar{\varphi}\bar{\varphi}$ is given by Eq. (51). Thus, the upper bound on the conversion cross section (33) can be rewritten as an upper bound on the quartic coupling:

[TABLE]

Note, however, that this constraint becomes less severe if the initial asymmetry is generated at larger values of $x_{0}>1$, see also discussion in Section III. We conclude that the exact bound on $\lambda_{\rm conv}$ is very sensitive to the moment $x_{0}$ and cannot be predicted purely within CP4 3HDM.

If the value of $\sigma_{\rm conv}$ is small compared to the upper bound (34) and $\delta_{0}$ is not vanishingly small, then $\delta$ approaches unity during the recovery stage, which modifies the freeze-out moment $x_{f}$ shifting it to smaller values. The relic density in this case is determined by the value of the absolute asymmetry $\Delta_{*}$ at the moment of the effective decoupling of the conversion processes, which also sets the relic density:

[TABLE]

This result is interesting for two reasons. First, similarly to the asymmetric case without conversion, we see that the relic density does not depend on the value of $\sigma_{\rm SM}$, given that it is not extremely small so that the symmetric freeze out would occur later than the asymmetric one. Second, from the model-building perspective, the presence of the conversion brings more parametric freedom, since the value of the relic density is not fixed by the initial conditions.

In summary, depending on the strength of the conversion process, we expect to observe two distinct regimes. For sufficiently weak conversion, with the bound (34) satisfied, one expects that the initial asymmetry, evolving through a depleted intermediate stage, will recover by the freeze-out moment $x_{f}$ and can eventually grow to values close to unity. The resulting relative contribution of antiparticles can be extremely small if the annihilation into SM is effective so that the parameter $a\Delta_{*}$ is large enough, see Eq. (26). For a stronger conversion process, the recovery is only partial. Even if the conversion process by itself is very weak at the moment of freeze-out $x_{f}$, it was stronger at earlier times and it shifted the position of the minimum $x_{*}$ to larger values. As a result, $\delta$ did not have enough time to grow to unity, and its final value can be rather small.

III Asymmetry evolution: numerical analysis

Though there are several available codes for numerical computation of the relic abundance, including micrOMEGAs Belanger:2014vza , DarkSUSY Bringmann:2018lay and MadDM Ambrogi:2018jqj , none of them is capable of handling asymmetric DM models with particle-antiparticle conversion processes. For that reason we have developed our own C++ code.

We use two different integration methods at different stages of the thermal evolution: the implicit Euler method for $x\sim 1\div 20$ and the classical Runge-Kutta method (RK4) for x\mathrel{\hbox to0.0pt{\lower 4.0pt\hbox{\hskip 1.0pt\sim}\hss}\raise 1.0pt\hbox{>}}20. While the latter method has a better precision ${\cal O}(t^{4})$ for a given integration step $t$, it fails to provide adequate solutions when the total density $Y_{t}$ is extremely close to the equilibrium density $Y_{\rm eq}$. The backward Euler method is the simplest implicit method, which manages to solve the system (14) in this stiff regime. The global truncation error for this method is ${\cal O}(t)$, so one has to use very small integration steps to get sufficient precision. In order to validate the solution obtained with the backward Euler method, we applied the implicit midpoint method, which has an accuracy ${\cal O}(t^{2})$. For $t\lesssim 10^{-3}$, the results of these two methods are consistent up to ${\cal O}(t)$.

We also compared our results in the symmetric case and in the absence of the conversion interaction to a simple freeze-out analytic approximation of the relic density gorbunov2011introduction

[TABLE]

and to the results of micrOMEGAs. As an example, for $m=250\,{\rm GeV}$ and $\sigma_{0}=9.8\cdot 10^{-8}\,{\rm GeV}^{-2}$ the formula above gives $\Omega h^{2}=0.22$. In the case of the Singlet DM model, where the values of the parameters correspond to $m$ and $\sigma_{0}$ given above, micrOMEGAs returns $\Omega h^{2}=0.16$, while our code gives $\Omega h^{2}=0.2$.

For our analysis of the asymmetry evolution we adopt the following default values of the model parameters: $m=250\,{\rm GeV}$, $\lambda_{346}=0.2$, $\lambda_{3}=0.9$, $\lambda_{\rm conv}=10^{-5}$. In the model we consider, which is presented in detail in Appendix A, the second and the third parameters determine the intensity of the interaction with the SM sector. The parameter values presented above correspond to $\langle\sigma v\rangle_{\!\,{\rm SM}}\approx\sigma_{\rm SM}\sim 10^{-7}\,{\rm GeV}^{-2}$. Given that $\langle\sigma v\rangle_{\!\,{\rm SM}}$ is much larger than $\langle\sigma v\rangle_{\rm conv}$, the last parameter in the list affects the asymmetry evolution the most, so our main goal is to study its variations.

When implementing the evolution of asymmetric DM, we need to fix the initial conditions: that is, we need to select the starting point $x_{0}$ and select the initial total density $Y_{t0}$ and the initial asymmetry $\delta_{0}$. Our default choice is $x_{0}=1$, $\delta_{0}=0.1$ and $Y_{t0}=Y_{\rm eq}$. This initial choice brings us from the start sufficiently close to equilibrium for the total density evolution. As a stability check, we also performed runs in which $\delta_{0}=0.1$ but the initial total density was chosen to be $10Y_{\rm eq}$ and $100Y_{\rm eq}$, which are very out-of-equilibrium situations. In all the cases, we observe (see Fig. 2) that the total density approaches the equilibrium density very fast, within $\Delta x\ll 1$. This quick equilibration of $Y_{t}$ also produces the very sharp jump of $\delta$ to the vicinity of one; although the three curves at the lower plot of Fig. 2 looks as if starting from different values of $\delta$, the initial value in all of them is $\delta_{0}=0.1$. The small discontinuity of the curves visible in Fig. 2 at $x=1.25$ is the artifact of our code, in which the number of relativistic degrees of freedom $g$ changes discretely with $x$ and one of these jumps takes place at $T=200\,{\rm GeV}$. This is the only plot, where the minimal step between points in the interval $x=\left[1,2\right]$ is taken to be $10^{-3}$ to show in detail the transition from the initial state to the equilibrium.

Fig. 3 shows the main feature we discussed above: the $x$ evolution of the absolute asymmetry $\Delta$ and the relative asymmetry $\delta$ as well as the modifications induced by the conversion process $\varphi\varphi\leftrightarrow\bar{\varphi}\bar{\varphi}$. One sees that when $\lambda_{\rm conv}$ is sufficiently small (red curve), the asymmetry first drops but then quickly restores to unity, and the relic density freezes out at values of $x$ smaller than the typical symmetric DM value $x\approx 25$. For $\lambda_{\rm conv}\approx 10^{-5}$ (yellow curve), the minimum gets deeper, nevertheless the asymmetry fully recovers at later stages, around $x\approx 30$. This is the boundary value at which the recovery completes just in time when the total density freeze out takes place.

For larger values of the conversion coupling (green curve), depletion of the asymmetry is more intense and extends to slightly larger values of $x$. The recovery stage starts with a much lower value of $\delta$ and it does not reach unity by the moment, when the total density freezes out. There is some residual growth of $\delta$ in the flattening stage, but the final value of $\delta(\infty)\ll 1$. Taking $\lambda_{\rm conv}$ to be twice larger than the boundary value (magenta curve) induces conversion effects so strong that they irreparably wash out the asymmetry. Also note that the left part of Fig. 3 shows very nice correspondence between the numerical results and analytical predictions for the recovery stage, described by Eq. (19).

In Fig. 4, left, we show the role of the starting value of the asymmetry $\delta_{0}$ at $x_{0}=1$ for the fixed conversion strength $\lambda_{\rm conv}=10^{-5}$. One can see the neat proportionality of the curves, which confirms the results of the qualitative analysis: during each stage, the value of asymmetry $\delta$ drops or increases by certain factor, unless it approaches unity. This proportionality is driven essentially by the fact that $d\delta/dx$ is proportional to $\delta$ itself.

These three plots, Figs. 3 and 4, illustrate our main conclusion: in order for the asymmetry to get restored close to 100%, one needs to select not too small $\delta_{0}$ and not too strong $\lambda_{\rm conv}$. Both of these values are the free parameters of the model, and Fig. 4, right, demonstrates the boundary between the two regimes on their plane.

The exact boundary, of course, depends on the moment $x_{0}$ at which certain initial asymmetry $\delta_{0}$ is postulated. In all plots so far we assumed that $x_{0}=1$. If the asymmetry is generated at smaller temperatures (larger values of $x_{0}$), then the asymmetry can get fully restored even in the presence of more intense conversion process. This dependence is rather sharp and can be observed already in Fig. 4, left. If one starts at $x_{0}=1$ with $\delta_{0}=10^{-3}$ (magenta curve), then $|\lambda_{\rm conv}|=10^{-5}$ keeps the final $\delta$ at the few percent level. However, if one postulates the same $\delta_{0}=10^{-3}$ at slightly lower temperature, with $x_{0}=1.5$, then the final asymmetry easily recovers to nearly unity. This strong sensitivity to the initial conditions is unavoidable when the depletion and recovery factors at the intermediate stages display exponential behavior.

IV Indirect detection

In the classical symmetric DM scenario, the single process — DM annihilation into the SM fields — determines, albeit in different kinematics, both the relic abundance and the possible indirect detection (ID) signal, which may be observable today. In the asymmetric DM scenarios, this relation is looser. In particular, in the highly asymmetric case of CP4 3HDM, even if there remains a vanishingly small amount of $\bar{\varphi}$ after freeze-out, it can be regenerated later via $\varphi\varphi\to\bar{\varphi}\bar{\varphi}$ conversion, and only then annihilate to the SM fields. In this Section, we investigate whether this two-stage process can produce an observable ID signal.

Let us assume that the relic abundance of $\bar{\varphi}$ is much smaller than the abundance of $\varphi$. Then, in the cosmologically modern epoch, the evolution of the local density of $\bar{\varphi}$ is given by the following equation

[TABLE]

It is implicitly assumed that all the quantities in Eq. (37) including the velocity-averaged cross sections can be spatially dependent. The density of $\varphi$ can be considered as approximately constant in time. Indeed, even if we consider symmetric DM in very dense regions, such as DM cusps at the centers of galaxies, where the density can be up to $\rho\sim 10^{3}\,{\rm GeV}/\,{\rm cm}^{3}$ (see e.g. Fig. 1 in Cirelli:2010xx ), and take the maximal value of $\langle\sigma v\rangle_{\!\,{\rm SM}}\sim 10^{-20}\,{\rm cm}^{3}/\,{\rm s}$ allowed by the unitarity bound for $m\sim 100\,{\rm GeV}$, the characteristic timescale of annihilation $\tau=\left(n\langle\sigma v\rangle_{\!\,{\rm SM}}\right)^{-1}$ would be around $10^{3}$ Gyrs, which already exceeds the lifetime of the Universe by two orders of magnitude. Our aim is to understand whether this large and essentially constant amount of DM can lead to efficient conversion and observable ID signal.

Taking into account our assumption, we can easily solve Eq. (37) for $\bar{n}$:

[TABLE]

Here, $\bar{n}_{0}$ is the remaining density of antiparticles after the freeze out, which can be extremely low as it is exponentially sensitive to the cross section of annihilation into SM, Eq. (26). The time $t$ is of the order of the Universe’s lifetime, $t\approx 10$ Gyrs, and it is much smaller that $\tau$. This allows us to simplify the solution for $r\equiv\bar{n}/n$ to

[TABLE]

Notice that $\langle\sigma v\rangle_{\rm conv}$ is calculated here in the present cosmological epoch. Due to its proportionality to the relative velocity, it is roughly two to three orders of magnitude smaller than it was at the freeze-out time (the exact value would be sensitive to the details of the DM velocity distribution in galaxies).

The bracket here can have either sign, depending on $r_{0}$, the relative amount of $\bar{\varphi}$ after the freeze out, and on the ratio of the two cross sections. Consider first the case when the sign is negative, which can happen if $r_{0}$ is not negligibly small. It corresponds to the evolution scenario where the relative amount of $\bar{\varphi}$ was initially fixed at the level of Eq. (26) after the freeze-out process, but then, in the structure formation epoch, the annihilation gradually switched back on, leading to further depletion of $\bar{\varphi}$’s and producing ID signal. In this regime, the amount of $\bar{\varphi}$ does not change much during the lifetime of the Universe: $r(t)\approx r_{0}$, and the conversion process does not play any role in forming the ID signal. The annihilation rate per unit volume is given by

[TABLE]

We now compare this with the rate of annihilation in the conventional symmetric case, which we write as $\gamma_{\rm sym}=n^{2}\langle\sigma v\rangle_{\!\,{\rm SM},0}$, with $\langle\sigma v\rangle_{\!\,{\rm SM},0}\approx 10^{-26}\,{\rm cm}^{3}/\,{\rm s}$ being the commonly adopted value of the velocity-averaged cross section that yields the correct value of the DM relic abundance. Here we would like to stress again that the relic abundance in the asymmetric scenario does not necessarily depend on $\langle\sigma v\rangle_{\!\,{\rm SM}}$. Thus the value of $\langle\sigma v\rangle_{\!\,{\rm SM}}$ can be larger than $\langle\sigma v\rangle_{\!\,{\rm SM},0}$ mentioned above and still be compatible with the observed relic abundance. Neglecting the possible difference in the spatial distribution of the annihilation signal that could appear in the presence of conversion, we get a simple expression

[TABLE]

Can this ratio be larger than 1? Naively, one may think of the situation when $r_{0}$ is only moderately small, while the SM annihilation cross section is much higher in our case than in the symmetric DM scenario. Unfortunately, this situation is not possible. The value of $r_{0}$, which is the same as $r(\infty)$ in Eq. (26), does depend on the SM annihilation (parameter $a$) in a very essential way. By increasing $\langle\sigma v\rangle_{\!\,{\rm SM}}$ with respect to $\langle\sigma v\rangle_{\!\,{\rm SM},0}$, one would strongly suppress $r_{0}$, deteriorating the ID signal.

One can also consider a different regime, when $r_{0}$ is extremely small due to a rather large annihilation cross section. Then, the sign of the bracket in Eq. (39) is positive. This is a novel scenario, in which particles $\bar{\varphi}$, whose amount was vanishingly small after the freeze-out, are being regenerated in the present cosmological epoch. To consider the extreme case, we set $r_{0}\to 0$ and obtain

[TABLE]

In this regime, the ID signal is time-dependent; it is predicted to be stronger now than during the structure formation epoch. However, its absolute value is even lower than in the previous case. An order-of-magnitude estimate for the values mentioned above gives

[TABLE]

The bottom line is that, for a very broad range of velocity-averaged conversion cross sections, the signal from dense DM regions in the asymmetric case is suppressed by many orders of magnitude compared to the symmetric DM scenario. Due to the strong constraints on the efficiency of the conversion process, we simply cannot generate any sufficient amount of $\bar{\varphi}$ in the given time.

V Conclusions

Dark matter can be stabilized by a $CP$ symmetry — provided it is a symmetry of an order higher than 2 and it is preserved by the vacuum. CP4 3HDM, the three-Higgs-doublet model featuring $CP$ symmetry of order 4, is the simplest framework producing such DM candidates Ivanov:2015mwl . Motivated by this form of stabilization, which was not observed in beyond the SM models before, we investigate its consequences for the DM evolution and freeze-out process.

Due to the unavoidable mass degeneracy Haber:2018iwr , the DM sector contains two stable complex fields $\varphi$ and $\bar{\varphi}$, which are CP4 eigenstates but which are not particle and antiparticle of each other. It is a key consequence of $CP$ conservation of this model that $\varphi\bar{\varphi}$ can annihilate into SM fields, while $\varphi\varphi$ and $\bar{\varphi}\bar{\varphi}$ cannot. However, these two pairs can turn into each other, $\varphi\varphi\leftrightarrow\bar{\varphi}\bar{\varphi}$, and this conversion process can strongly affect the DM evolution. It offers a rare example of a DM model in which self-interaction in the dark sector can significantly affect the cosmological observables and, therefore, can be constrained.

These considerations led us to investigate the asymmetric DM regime, in which some asymmetry between $\varphi$ and $\bar{\varphi}$ is generated at $T\sim m$ by a mechanism external to CP4 3HDM. The presence of conversion process makes this model clearly different from the more conventional asymmetric DM models Graesser:2011wi ; Iminniyaz:2011yp ; Zurek:2013wia ; Petraki:2013wwa . The questions we wanted to answer included: does the asymmetry survive in the presence of conversion? Is the final asymmetry large or small? What is the role of the initial conditions and of the physical processes? In the strongly asymmetric case, can this model produce an observable indirect detection signal through the two-stage conversion-then-annihilation process?

To answer them, we presented a detailed analytical and numerical analysis of the system of coupled Boltzmann equations, without and with conversion. For the numerical analysis, we wrote our own code since the existing packages are incapable of dealing with this model. We found that the relative asymmetry $\delta$ first sharply drops from its initial value, driven by the conversion process. When conversion effectively switches off, the relative asymmetry enters the recovery stage, and finally freezes out, either at some small value $\delta(\infty)\ll 1$, implying $n(\bar{\varphi})\approx n(\varphi)$, or very close to 1, which means $r=n(\bar{\varphi})/n(\varphi)\ll 1$. We established the boundary between the two regimes and showed how it depends on the initial conditions and on the intensity of the conversion process. In particular, we established that, in order for the asymmetry to recover to sizable values, the conversion cross section must be several orders of magnitude smaller than the annihilation one.

We also checked the intensity of the ID signal in both regimes and found that the expected signal is much smaller than in the symmetric DM scenario. Although it does not offer a novel observational signature of this model, it may help to resurrect scenarios in which one would normally expect a very strong ID signal in the usual symmetric DM case.

We stress that we did not aim to perform a systematic scan of the parameter space of CP4 3HDM. Nor do we claim that this model is superior to others in resolving some long-standing DM-related issues, such as small-scale problems. Rather, we wanted to investigate the phenomenon itself. We believe that qualitative features of the thermal DM history found in this model can arise in other models with the conversion processes competing with annihilation, which deserves further investigation.

I.P.I. acknowledges funding from the the Portuguese Fundação para a Ciência e a Tecnologia (FCT) through the FCT Investigator contract IF/00989/2014/CP1214/CT0004 under the IF2014 Programme, and through the contracts UID/FIS/00777/2013, UID/FIS/00777/2019, CERN/FIS-PAR/0004/2017 and PTDC/FIS-PAR/29436/2017, which are partially funded through POCTI (FEDER), COMPETE, QREN, and the EU. I.P.I. also acknowledges support from National Science Center, Poland, via the project Harmonia (UMO-2015/18/M/ST2/00518). I.P.I. thanks University of Liege for hospitality and stimulating atmosphere. M.L. thanks A. Bhattacharya for useful advice on C++. The work of M.L. is supported by a FRIA grant (F.N.R.S.).

Appendix A DM states and interactions in CP4 3HDM

The two lightest neutral fields $h$ and $a$ in Eq. (3) are the DM candidates protected by CP4 symmetry. In principle, all interactions can be written in terms of these real fields. However, they are not CP4-eigenstates, as they transform under CP4 as $h\xrightarrow{CP}-a$ and $a\xrightarrow{CP}h$. Since CP4 is conserved, we combine them into CP4 eigenstates Ivanov:2015mwl ; Aranda:2016qmp :

[TABLE]

We stress that this peculiar form of CP transformation does not involve any redefinition of the CP symmetry per se; it is the same CP4 as in Eq. (2), just applied to the physical states after minimization Aranda:2016qmp .

Much care must be taken when interpreting particle/antiparticle assignments of the corresponding one-particle states. Express the complex field $\varphi$ via creation and annihilation operators,

[TABLE]

where $px\equiv Et-\vec{p}\vec{x}$ and $\tilde{dp}\equiv d^{3}p/[2E(2\pi)^{3}]$. The standard assignment of the $C$-transformation, which would map a complex field to its conjugate, would interchange operators $a$ and $b$ and would make the one-particle states $a^{\dagger}|0\rangle$ and $b^{\dagger}|0\rangle$ antiparticles of each other. In our case, the definition of the $CP$-symmetry via Eq. (44) forces us to accept that the one-particle state $a^{\dagger}|0\rangle$, which we denote as $\varphi$ in the main text, is its own antiparticle. The other one-particle state with the same mass, $b^{\dagger}|0\rangle$, which we denote as $\bar{\varphi}$, is a different state and is also its own antiparticle, see discussions in Aranda:2016qmp ; Haber:2018iwr . We stress again that we abused the language in the main text of this paper and called $\bar{\varphi}$ “antiparticles”. This was done only to facilitate wording; in reality, $\varphi$ and $\bar{\varphi}$ are not antiparticles of each other, but are two distinct mass-degenerate particles.

The full Lagrangian can be rewritten in terms of these complex fields rather than $h$ and $a$ Kopke:2018vyw ; Haber:2018iwr . In this paper, we do not perform a detailed scan of the parameter space of the model but rather focus on a particular phenomenon. Thus, we find it sufficient to consider the DM mass above 100 GeV and to consider, for annihilation, only the channels $\varphi\bar{\varphi}\to WW$, $ZZ$. Their diagrams include the pointlike interactions and the $s$-channel annihilation via the SM-like Higgs boson coming from

[TABLE]

where $\bar{g}^{2}\equiv g^{2}+g^{\prime 2}$, with $g$, $g^{\prime}$ being the gauge coupling constants, as well as the $t$ and $u$-channel annihilation diagrams which involve triple interactions linking $\varphi$’s with the gauge bosons and the heavier inert scalars, either neutral or charged. In numerical calculations, we used all these diagrams. To show an estimate of the cross section we expect at small incoming velocities and large $m\gg m_{W}$, one can neglect the $t$ and $u$ channel diagrams and express the cross section of the DM annihilation into the SM fields as

[TABLE]

Here, the first two contributions in the brackets come from the transverse polarizations of the $W$ and $Z$ pairs; they combine to $\approx 0.08$. The last contribution is the large-$m$ estimate of the cross section of the longitudinal gauge boson pairs.555For a pedagogical exposition, with the example of the inert doublet model, of how $\lambda_{3}$ arises in the annihilation amplitude, see Section 3.2 of Ref. Hambye:2009pw Thus, for $\lambda_{3}=0.9$ and $m=250$ GeV used in numerical calculations, we get

[TABLE]

Of course, in the funnel region, $m\approx m_{h_{125}}/2$, the annihilation cross section will be resonantly enhanced, while at even smaller masses it will be strongly suppressed.

For the conversion process, we first write the self-interaction potential for $\varphi$:

[TABLE]

The $\varphi^{4}$ interaction terms lead to the following amplitude of conversion process $\varphi\varphi\to\bar{\varphi}\bar{\varphi}$:

[TABLE]

whose cross section is

[TABLE]