Phases of Dark Matter from Inverse Decays

Abstract

Inverse decays are an interesting avenue for producing dark matter in the early universe. We study in detail various phases of dark matter parameter space where inverse decays control its abundance, expanding on our work of INDY dark matter and going beyond. The role of initial conditions and the impact of departure from kinetic equilibrium are investigated as well. We show how these inverse decay phases can arise in theories of a kinetically mixed dark photon and dark Higgs, with promising prospects for detection at upcoming experiments.

I Introduction

The identity of dark matter (DM) remains one of the most pressing puzzles of modern day physics. Recent years have seen a surge of new ideas suggesting various processes in the early universe that could control the DM relic abundance, including SIMPs Hochberg et al. (2014, 2015), ELDERs Kuflik et al. (2016), forbidden channels Griest and Seckel (1991); D’Agnolo and Ruderman (2015), coscattering D’Agnolo et al. (2017), zombies Kramer et al. (2021) and more, where in some cases decay interactions play a role in the evolution over time of the DM abundance Dror et al. (2016, 2018); Kopp et al. (2016); D’Agnolo et al. (2018); Fitzpatrick et al. (2022); Kim and Kuflik (2019); Asadi et al. (2021a, b); Dror et al. (2016, 2018); Berlin et al. (2016); Morrissey et al. (2009); Cohen et al. (2010); Bandyopadhyay et al. (2011); Farina et al. (2016); Kim and Kuflik (2019); Hochberg et al. (2019); Feng et al. (2003); Kaplinghat (2005); Moroi and Randall (2000); Acharya et al. (2009); Hall et al. (2010); Berlin et al. (2016); Morrissey et al. (2009); Azatov et al. (2021); Garny et al. (2017); Yaguna (2008); Garny et al. (2019). Among these, an exciting possibility is for inverse decays of the DM to control its relic abundance. Ref. Frumkin et al. (2023a) proposed a new thermal DM candidate whose abundance is determined by the freezeout of inverse decays, dubbed ‘INDY’ DM. In this work, we elaborate on the INDY mechanism and study in detail the various phases that arise when considering such inverse decays.

This paper is organized as follows. We begin in Section II by setting up the basic mechanism and in Section III we identify the relevant theory phases. In Section IV we scrutinize the effects of departure from kinetic equilibrium. Section V presents a model that realizes these phases, along with its phenomenological prospects. We conclude in Section VI.

Figure 1: Phases of DM. Annihilations versus decays phase diagram, for various mass splittings and DM masses, and different initial conditions. The solid lines in the right (where $\alpha_{\rm decay}>10^{-24}$) are the same as in Ref. Frumkin et al. (2023a), with initial condition of chemical equilibrium at $x=1$. The two phases are evident: The first is the horizontal branch, where the relic abundance of the DM is set through the annihilations of another dark particle, with which it is in chemical equilibrium via decays and inverse decays, denoted by ‘coannihilation’. The second is the vertical branch, where the relic abundance of the DM is set by the freezeout of the inverse decay of the DM particle. We dub this case as ‘INDY’ DM. The lighter-colored solid lines in the right represent the ‘Freeze-in-Freeze-out’ (FIFO) case, where the initial conditions taken are $Y_{\chi}=0$ at $x=0$, and the $\chi$ abundance grows and decreases. The FIFO case coincide with INDY for most of the parameters, or otherwise very close to it. In the left side of the figure ($\alpha_{\rm decay}>10^{-24}$) the lines represent the ‘Freeze-in’ (FI) case, where the initial conditions taken are $Y_{\chi}=0$ at $x=0$, and the $\chi$ abundance grows with time, as well as the ‘FO & Decay’ phase, where $\psi$ freezes out and later decays to the DM.

II The Basic System

We begin by setting up the system of interest. Consider the decay and inverse decay processes $\psi\longleftrightarrow\chi+\phi$ of a DM candidate $\chi$ and an unstable particle $\psi$, along with the self-annihilation process $\psi\psi\longleftrightarrow\tilde{\phi}\tilde{\phi}\$, with $\phi$ and $\tilde{\phi}$ indicating particles in equilibrium with the Standard Model (SM) bath. The relative strength between these two processes will determine which reaction controls the relic abundance of DM $\chi$ and will lead to various phases of this simple system.

The Boltzmann equations (BEs) for this system, assuming detailed balance, can be written as Frumkin et al. (2023a)

$$\frac{dY_{\chi}}{dx} = -\frac{\langle\Gamma\rangle}{xH}\left(\frac{Y_{\psi}^{\rm eq}}{Y_$$(1)

$$\frac{dY_{\psi}}{dx} = -\frac{dY_{\chi}}{dx}-\frac{s\langle\sigma v\rangle}{xH}\left(Y_{$$

where $x=m_{\chi}/T$, $Y=n/s$ is the yield, $s$ is the entropy density, $H$ is the Hubble rate, $\langle\Gamma\rangle$ represents the thermally averaged (time dilation included) decay rate of $\psi\to\chi\phi$, and $\langle\sigma v\rangle$ represents the thermally averaged cross section for the annihilation $\psi\psi\to\tilde{\phi}\tilde{\phi}$. (The changes in the effective number of degrees of freedom for radiation density $g_{*}$ and for entropy $g_{*s}$ are taken into account by dividing $H$ by a factor of $\left(1+\frac{1}{3}\frac{d{\rm ln}g_{*s}(T)}{d\ln T}\right)$; values are taken from Ref. Hindmarsh and Philipsen (2005).)

We parameterize the decay width and the thermally averaged cross-section of the $\psi$ self-annihilations by

$$\Gamma\equiv\alpha_{\rm decay}m_{\psi}\bar{\beta}$$(2)

and

$$\langle\sigma v\rangle\equiv\frac{\alpha_{\rm ann}^{2}}{m_{\psi}^{2}}\,$$(3)

respectively, where $\alpha_{\rm decay}$ and $\alpha_{\rm ann}$ are dimensionless couplings constants we use to parameterize the system and $\bar{\beta}=\sqrt{1-2\frac{\left(m_{\phi}^{2}+m_{\chi}^{2}\right)}{m_{\psi}^{2$ is a phase space factor. The thermally averaged decay rate obtained by using Maxwell-Boltzmann statistics is given by

$$\left\langle\Gamma\right\rangle=\Gamma\frac{K_{1}\left(\frac{m_{\psi}}{T}$$(4)

with $K_{1,2}$ the modified Bessel functions of the second kind.

In what follows we use

$$\Delta\equiv\frac{m_{\psi}-m_{\chi}}{m_{\chi}}\quad\text{and}\quad r\equiv$$(5)

to denote the masses $m_{\psi}$ and $m_{\phi}$ compared to the mass of the DM $m_{\chi}$. We further take the number of degrees of freedom of the dark particles to be $g_{\chi}=g_{\psi}=4$ for concreteness. The BEs (1) depend on $m_{\phi}$ only through the decay rate $\Gamma$; we thus take $r=0$ in upcoming sections unless specified otherwise, and note that the massive $\phi$ case (i.e.$r\neq 0$) can be inferred by augmenting the coupling $\alpha_{\rm decay}$ of Eq. (2) by the ratio of the relative allowed phase spaces.

The BEs (1) admit several distinct phases that reproduce the observed DM relic abundance, depending on masses, couplings, and initial conditions. These phases are illustrated in Fig. 1, showing the values of couplings $\alpha_{\rm decay}$ and $\alpha_{\rm ann}$ that reproduce the observed abundance of DM, $Y_{\rm obs}$, for various cases. Throughout this paper, we parameterize the observed DM abundance as $Y_{\rm obs}=c\;T_{\rm eq}/m_{\chi}$, with $T_{\rm eq}=0.8$ eV the temperature at matter-radiation equality and $c=0.54$Hochberg et al. (2014).

III Phases

We now delve into each of the phases that appear in Fig. 1.

Figure 2: BE evolution. Examples of solutions to the BEs  (1) that reproduce the observed DM relic abundance for different parameters. In each panel, we present the calculated (solid) and equilibrium (dashed) abundances for $\chi$ (blue) and $\psi$ (orange) as a function of $x$, along with the observed relic abundance (green). Both panels represent the vertical branch in which the freeze out of inverse decay determines the relic abundance. The left panel represents coannihilations via decay. The middle panel corresponds to INDY DM that is out of chemical equilibrium, while the right panel corresponds to an INDY that is in chemical equilibrium.

III.1 Coannihilating via decays

First, we consider the phase obtained for large $\alpha_{\rm decay}$, which is demonstrated by the horizontal curves of constant $\alpha_{\rm ann}$ towards the right side of Fig. 1, where $\alpha_{\rm decay}>10^{-24}$. Here, the relic abundance of $\chi$ is being set by the annihilations of $\psi$ into other particles. The rapid decay and inverse decay rates, due to the large coupling $\alpha_{\rm decay}$, keep the particles $\chi$ and $\psi$ in chemical equilibrium (CE) with each other, balancing their number density ratios to be as in equilibrium $n_{\chi}/n_{\psi}=n_{\chi}^{\rm eq}/n_{\psi}^{\rm eq}$ Griest and Seckel (1991). This enables one to write the system as a single BE for the sum of the number densities,

$$\frac{dn}{dt}+3Hn=-\frac{\langle\sigma v\rangle}{\left(1+\frac{n_{\chi}^{\rm eq$$(6)

with $n\equiv n_{\psi}+n_{\chi}$. Eq. (6) is general for any strongly coupled two-particle system and is not unique for inverse decays. Indeed, a similar equation is found for the well-studied case of coannhilations Griest and Seckel (1991); the difference is that here the CE between the two particles is maintained via decays and inverse decays, rather than annihilations.

The instantaneous freezeout approximation provides an estimated solution for the requisite cross section explaining the observed DM abundance, or said another way, provides a relationship between the mass $m_{\chi}$ and the annihilation coupling $\alpha_{\rm ann}$:

$$m_{\chi}\simeq\alpha_{\rm ann}^{\frac{2}{2+\Delta}}\left(m_{\rm pl}T_{\rm eq}^$$(7)

where ${m_{\rm pl}=2.4\times 10^{18}}$ GeV is the reduced Planck mass Griest and Seckel (1991); Frumkin et al. (2023b). The left panel of Fig. 2 shows an example of the evolution of $Y_{\chi}$ and $Y_{\psi}$ over time along this phase. Here both particles are in equilibrium before freezeout occurs, since $\psi$ is coupled to the SM via annihilation and $\chi$ via strong coupling to $\psi$, and both particles depart from equilibrium simultaneously when the effective annihilation reaction rate is comparable to the Hubble rate, $n\langle\sigma v\rangle/\left(1+{n_{\chi}^{\rm eq}}/{n_{\psi}^{\rm eq}}\right)$.

The power law relation of Eq. (7) agrees well with numerical results, as is shown in the left panel of Fig. 3. The right panel of Fig. 3 shows the mass splitting $\Delta$ versus the DM mass $m_{\chi}$ at fixed values of $\alpha_{\rm ann}$ that reproduce the observed DM yield. $\alpha_{\rm ann}$ is generally constrained by unitarity arguments, which in turn places an upper limit on the possible DM masses along this phase.

III.2 INDY DM

We now move to the phase obtained for large $\alpha_{\rm ann}$ values, demonstrated by the vertical curves of constant $\alpha_{\rm decay}>10^{-24}$ towards the right side of Fig. 1. In this phase, first proposed in Ref. Frumkin et al. (2023a), the relic abundance of $\chi$ is set by the freeze out of inverse decays. We dub this ‘INverse DecaY’ (INDY) dark matter.

For large values of $\alpha_{\rm ann}$ (compared to the values obtained in the coannihilating-via-decays phase), $\psi$ preserves CE throughout the relevant times, and the BEs for the yield can be approximated by a simple linear differential equation,

$$\frac{dY_{\chi}}{dx}+\mathbf{a}(x)Y_{\chi}=\mathbf{b}(x)\,,$$(8)

with

$$\mathbf{a}(x)\equiv\frac{\langle\Gamma\rangle}{xH}\frac{Y_{\psi}^{\rm eq}}{Y_{$$(9)

denoting the (normalized) inverse decay and decay rates for $Y_{\chi}$, respectively.

Figure 3: Coannihilating via decays.Left: The annihilation coupling that reproduces the observed DM abundance along the horizontal phase as a function of DM mass, for various mass splittings. We show the numerical solution (solid), the instantaneous freezeout Kolb and Turner (1990) (dotted) and Griest-Seckel Griest and Seckel (1991) (dashed) approximations. Right: The mass splitting that reproduces the observed DM abundance along the coannihilation-via-decay phase as a function of DM mass, for fixed annihilation couplings. It is evident that increasing the mass or the mass splitting increases the required annihilation strength.

Eq. (8) can be integrated into the form

$$Y_{\chi}(x)=e^{-\mathbf{A}(x)}Y_{\chi}\left(1\right)+\intop_{1}^{x}e^{\mathbf{$$(10)

where the initial condition used throughout this section is an equilirubium yield at $T=m_{\chi}$, namely $Y_{\chi}\left(1\right)=Y_{\chi}^{\rm eq}\left(1\right)$ and

$$\mathbf{A}(x)\equiv\intop_{1}^{x}\mathbf{a}(\xi)d\xi\,.$$(11)

The relic abundance is obtained by taking the limit ${x\rightarrow\infty}$ in Eq. (10), producing

$$Y_{\chi,\infty}=e^{-\mathbf{A}_{\infty}}Y_{\chi}\left(1\right)+\intop_{1}^{$$(12)

where $\mathbf{A}_{\infty}\equiv\mathbf{A}(\infty)$.

Eq. (12) can be understood intuitively. For a short interval $dx$, the probability of a $\chi$ particle to undergo inverse decay is $\mathbf{a}(x)dx$. Therefore, the probability for particles at $x=1$ not to inverse decay at all is given by $e^{-\intop_{1}^{\infty}\mathbf{a}\left(\xi\right)d\xi}$. Decays of $\psi$ are another source of $\chi$ particles. For a short interval $d\eta$ at the time $\eta$, the number of $\chi$ particles created by $\psi$ decays is $\mathbf{b}(\eta)d\eta$, and the probability of each not to inverse decay again is given by $e^{-\intop_{\eta}^{\infty}\mathbf{a}\left(\xi\right)d\xi}$. The total contribution of the decays is given by summing over all relevant times from $\eta=1$ to $\eta=\infty$, which explains Eq. (12).

The two sources of $\chi$ particles—from early times $x<1$ (through initial conditions) and from $\psi$ decays at later times $x>1$—yield two types of INDY dark matter. When the contribution from the former dominates the relic abundance, namely the $Y_{\chi}\left(1\right)$ term in Eq. (12) is dominant, the INDY DM candidate is out of CE. When the contributions from the subsequent $\psi$ decays are most relevant, the second $\mathbf{b}(x)$ term in Eq. (12) dominates the relic, and the INDY DM candidate follows CE. For a given DM mass and splitting, we determine numerically what kind of INDY DM arises by examining the ratio $w$,

$$w\equiv\frac{\intop_{1}^{\infty}e^{\mathbf{A}(\eta)}\mathbf{b}(\eta)d\eta}{Y_{$$(13)

and establishing whether $w>1$ (in CE) or $w<1$ (out of CE).

The left panel of Fig. 4 shows $\alpha_{\rm decay}$ as a function of $m_{\chi}$ for various mass splittings $\Delta$, where the observed abundance is obtained and DM is an INDY. The right panel of Fig. 4 depicts $w$ along the DM solutions, indicating that for small mass splittings (e.g.$\Delta=0.05$) the initial $\chi$ particles dominate the relic abundance, while for large mass splittings (e.g.$\Delta=0.25$) the later $\psi$ decays are more significant contributors to the $\chi$ abundance.

We can understand the slope and shape of the curves in the left panel of Fig. 4 of the INDY DM solutions by examining when the INDYs are either in or out of chemical equilibrium, as we now discuss.

Figure 4: INDY DM.Left: The decay coupling that reproduces the observed DM relic abundance along the vertical branches (for $\alpha_{\rm decay}>10^{-24}$) as a function of DM mass, for various mass splittings. We show the numerical solution to the Boltzmann equations; the dashed region of the curves indicates where the effective coupling $\alpha_{\rm ann}$ that we use to parameterize the $\psi$ annihilation cross section exceeds $100$ along the horizontal phase. Right: The ratio $w$ in Eq. (13) as a function of DM mass, for fixed mass splitting, along the solutions of INDY DM shown in the left panel. The dashed gray curve indicates $w=1$, where $w>1$ (or $w<1$) indicates where $\psi$ decays (or initial conditions) dominate the $\chi$ abundance.

III.2.1 INDYs out of chemical equilibrium ($w<1$)

Here $\psi$ decays can be neglected, and the relic abundance of $\chi$ is determined by the $Y_{\chi}\left(1\right)$ term in Eq. (12). $Y_{\chi}$ departs from chemical equilibrium early on, and the reaction rate is maximal at $x_{*}=\Delta^{-1}$. This can be seen by approximating $\mathbf{a}(x)$ in the non-relativistic (NR) limit as

$$\mathbf{a}(x)\approx\frac{3}{\pi}\sqrt{\frac{10}{g_{*}}}\frac{g_{\psi}}{g_{$$(14)

which is small at $x=1$, and so the reaction rate is not fast enough to bring $\chi$ toward equilibrium. A visual realization of such behavior is shown in the middle panel of Fig. 2.

The approximated relic abundance in this case is given by $Y_{\chi,\infty}\approx e^{-\mathbf{A}_{\infty}}Y_{\chi}^{\rm eq}\left(1\right)$, where $\mathbf{A}_{\infty}$ is obtained by integrating Eq. (14),

$$\mathbf{A}_{\infty}\approx\frac{3}{\pi}\sqrt{\frac{10}{g_{*}}}\frac{g_{\psi}}{$$(15)

Since the logarithm of the ratio between the initial and observed abundance $\mathbf{A}_{\infty}\approx\ln\left(\frac{Y_{\chi}^{\rm eq}\left(1\right)}{Y_{$ is of order unity, using Eq. (15) we find that the decay coupling which reproduces the observed abundance scales as

$$\alpha_{\rm decay}\sim\frac{m_{\chi}}{m_{\rm pl}}\,.$$(16)

This linear dependence between the $\alpha_{\rm decay}$ and $m_{\chi}$ explains the linear shape and slope of the small splitting $\Delta=0.05,\;0.1$ curves in the left panel of Fig. 4. The dependence on $\Delta$ comes only through the unwritten prefactor in Eq. (16) and not through the exponent of a mass scale as in Eq. (7). Note that since the relation in Eq. (16) is independent of $T_{\rm eq}$, the couplings that reproduce the DM relic abundance for this case are much smaller than in other mechanisms (compare e.g. Eq. (7)).

A comment is in order regarding initial conditions. Annihilation mechanisms (e.g. WIMP, SIMP and co-annihilation) are generally independent of initial conditions, since their reaction rate decreases with time up until freezeout occurs. However, for inverse decays initial conditions may play a role, since the reaction rate does not monotonically decrease. The potential initial condition dependence in this case is also manifested by the explicit dependence on $Y_{\chi}\left(1\right)$ in Eq. (12) for INDY dark matter.

The left panel of Fig. 5 demonstrates the above discussion. The colored curves show the relic abundance for a fixed DM mass $m_{\chi}=1$ GeV and splitting $\Delta=0.05$, corresponding to the solid orange curve in Fig. 1, as we vary the yield at $x=1$ from its equilibrium value. The larger values of $\alpha_{\rm decay}$ correspond to the coannihilating-via-decays phase, where the annihilation of $\psi$’s controls the $\chi$ abundance, and as expected is insensitive to the initial conditions. In contrast, for the smaller values of $\alpha_{\rm decay}$ that correspond to INDY DM, we see that the initial conditions can potentially impact the DM relic abundance.

The right panel of Fig. 5 shows the impact that the change in initial conditions has on the value of the coupling $\alpha_{\rm decay}$ along INDY DM. As is evident, while the relic abundance might vary as the initial conditions are changed, the change to the coupling required to reproduce the observed relic abundance is minimal. The reason can be traced to the exponential dependence on $\alpha_{\rm decay}$ when the INDY is out of chemical equilibrium, as discussed earlier: large changes in the prefactor can be compensated by small relative changes in the exponent, namely small changes to the coupling (see Eq. (12)).

We learn that the impact of initial conditions on the coupling $\alpha_{\rm decay}$ for INDY DM is at worst logarithmic in the ratio of the initial abundance to the equilibrium one.

Figure 5: INDY DM.Left: The relic abundance obtained versus the couplings presented in Fig. 1 for $m_{\chi}=1\ {\rm GeV}$ and $\Delta=0.05$, with different initial condition values. Right: The ratio of the decay couplings obtained for various ratios of initial conditions to equilibrium initial conditions, as a function of DM mass.

III.2.2 INDYs in chemical equilbrium ($w>1$)

For large values of $m_{\chi}$ and $\Delta$, the coupling $\alpha_{\rm decay}$ is large enough such that the inverse decay rate $\mathbf{a}(x)$ keeps the abundance $Y_{\chi}$ close to its equilibrium value at $x\sim 1$. This can be seen for example in the right panel of Fig. 2. Thus, the vast majority of the initial $\chi$ particles at $x=1$ are removed, and the relic abundance is controlled by the $\psi$ decay term in Eq. (12).

In the NR limit,

$$\mathbf{b}(x)=\frac{3^{3}5^{3/2}}{2^{2}\pi^{9/2}}\frac{g_{\psi}\bar{\beta}$$(17)

By using the saddle point approximation, as detailed in Appendix A, one can show that the relic abundance of DM scales with DM mass and decay coupling (up to logarithmic corrections) as

$$Y_{\chi,\infty}\sim\left(\frac{\alpha_{\rm decay}m_{\rm pl}}{m_{\chi}}\right)^$$(18)

Matching the relic abundance in Eq. (18) with the observed value yields the relation

$$\alpha_{\rm decay}\sim\frac{m_{\chi}^{1+\Delta}}{T_{\rm eq}^{\Delta}m_{\rm pl}$$(19)

This explains the linear shape and slope of the larger mass splitting $\Delta\geq 0.1$ curves in Fig. 4. In INDY dark matter, even for larger DM masses one obtains small decay couplings, significantly smaller than couplings for the WIMP.

III.3 FI and FIFO

As shown in Ref. Hall et al. (2010), there is theoretical motivation to consider DM with a negligibly small abundance at very early times. This case, which corresponds to initial conditions of $Y_{\chi}\left(0\right)=0$, leads to two more possible classes of solutions of the BEs (1), as indicated by the vertical curves in Fig. 1 which are not the INDY phase.

The left panel of Fig. 6 presents the relic abundance obtained for different $\alpha_{\rm decay}$, in the limit of large annihilation values (i.e$\psi$ is in CE), for initial conditions of $Y_{\chi}\left(0\right)=0$ (blue). As a reference, the INDY case, $Y_{\chi}\left(1\right)=Y_{\chi}^{\rm eq}\left(1\right)$, is plotted as well (orange). As can be seen, there are two possibilities of $\alpha_{\rm decay}$ values which reproduce the relic abundance, beyond the INDY solution. The first, indicated by the crossing of the black and blue curves at small $\alpha_{\rm decay}$, is a freeze-in (FI) mechanism Hall et al. (2010). The second, indicated by the crossing of the black and blue curves at larger values of $\alpha_{\rm decay}$, presents the case of the DM abundance growing at first and then decreasing—for some masses never tracking the equilibrium abundance. This is a freeze-in and freeze-out (FIFO) scenario. Examples of the evolution over time of the DM abundance for the cases of FI and FIFO DM are presented in Ref. Frumkin et al. (2023a).

Figure 6: FIFO, FI and INDY.Left: The relic abundance produced versus the decay coupling for different initial conditions for $Y_{\psi}=Y_{\psi,{\rm eq}}$. The intersections between the calculated curves (blue and orange) and the observed abundance line (black) provide the relevant couplings for the FI, FIFO and INDY cases. Right: The ratio between the decay couplings, which reproduce the observed abundance for FIFO and INDY cases, as a function of DM mass, for different mass splittings.

For FI or FIFO to occur for a given DM mass and mass splitting, there must exist a value of $\alpha_{\rm decay}$ which is large enough to create a sufficient number of $\chi$ particles (through $\psi$ decays), but small enough so that inverse decays do not eliminate all DM particles and the observed abundance remains. To check when these conditions can be met in the high $\alpha_{\rm ann}$ limit, one can examine the maximal value of $Y_{\chi,\infty}$ as a function of $\alpha_{\rm decay}$, with initial conditions of $Y_{\chi}\left(0\right)=0$. By taking the relic abundance obtained in this case,

$$Y_{\chi,\infty}=\intop_{0}^{\infty}e^{-\intop_{\eta}^{\infty}\mathbf{a}\left($$(20)

matching its derivative with respect to $\alpha_{\rm decay}$ to zero,

$$\intop_{0}^{\infty}\left[e^{-\intop_{\eta}^{\infty}\mathbf{a}\left(\xi\right)d$$(21)

solving for $\alpha_{\rm decay}$ and comparing the value of the relic abundance obtained with this $\alpha_{\rm decay}$ to $Y_{\rm obs}$, it is possible to determine if FI or FIFO can occur for each choice of mass and mass splitting. Visually, this condition is equivalent to requiring that the blue and black curves in Fig. 6 intersect. Numerically we find that the condition for either FI or FIFO is not fulfilled when $\Delta<0.1$ and $m_{\chi}\sim{\rm keV}$.

A comparison between the decay couplings yielding the FIFO case $\alpha_{\rm FIFO}$ and those yielding INDY DM $\alpha_{\rm INDY}$ is shown in the right panel of Fig. 6. As expected, the ratio is close to unity, in agreement with the discussion of the previous section that the coupling for freeze out has only a small dependence on initial conditions. The dips in the curves at ${\cal O}({\rm MeV})$ and $100~{}{\rm MeV}$ to GeV DM masses correspond to variations in $g_{*}$, dominantly from electron-positron annihilation and the QCD phase transition, respectively. On the other hand, the decays couplings for the FI case are almost independent of the mass scale, as in Ref. Hall et al. (2010). This is because for small couplings, as occurs in FI, the inverse decays can be neglected, and the relic abundance in Eq. (20) is approximated to

$$Y_{\chi,\infty}\approx\intop_{0}^{\infty}\mathbf{b}(\eta)d\eta\,.$$(22)

Since both $Y_{\rm obs}$ and $\mathbf{b}(\eta)$ are explicitly inversely proportional to $m_{\chi}$, the coupling values depend on the mass scale only via $g_{*}$ and $g_{*s}$.

III.4 Freezeout and Decay

The last phase presented in Fig. 1, which is manifested by the horizontal curves towards the left, is obtained by extremely low $\alpha_{\rm decay}$ values, and initial conditions of $Y_{\chi}\left(0\right)=0$. In this case, the $\psi$ annihilation is entirely identical to the WIMP, and the abundance of $\chi$ before and during the $\psi$ freezeout is negligible. Once $\psi$ freezes out, it decays into $\chi$, providing a similar relic abundance to the corresponding WIMP case, up to a factor of $\left(1+\Delta\right)^{-1}$.

Note that the monotonic increase of $\alpha_{\rm decay}$ with $\alpha_{\rm ann}$ for the curves in the left side of Fig. 1 is easily understood. For higher $\alpha_{\rm ann}$ values, fewer $\psi$ particles remain after their freezeout. To compensate, more $\chi$ particles must be created before $\psi$’s freezeout, and thus $\alpha_{\rm decay}$ must be larger. In this sense, this freezeout-and-decay phase is a continuation of the FI phase, with $\chi$ particles generated in different regions of the $\psi$ evolution.

IV Departure from Kinetic Equilibrium

In this section, we discuss non-kinetic equilibrium (NKE) corrections to INDY DM. Thus far, all the results shown and discussed were obtained by solving the integrated BEs (1) for various parameters and initial conditions. Using these equations, we assumed that $\chi$ is in kinetic equilibrium (KE) with the SM. While this assumption is valid when the decay and inverse decay rates are relatively large (in the coannihilation-via-decays phase), it is not necessarily valid for INDY DM, where the decay coupling is small and the corresponding rates are slow. (For a discussion of the impact of non-kinetic equilibrium, see e.g Refs. Binder et al. (2017); D’Agnolo et al. (2017); Garny et al. (2017).)

In this section, we relax the assumption of KE after $x=1$ and study its impact on INDY DM and the obtained couplings, in the limit of $Y_{\psi}=Y_{\psi}^{\rm eq}$ ($\alpha_{\rm ann}\rightarrow\infty$). We show that the solution of the NKE case provides small corrections to the values of $\alpha_{\rm decay}$ in the relevant parameter space, rendering the previously presented results valid. We further provide mathematical insight for these results. A reader interested directly in a QFT model realizing the phases discussed thus far can advance directly to the next section.

IV.1 Non-integrated Boltzmann Equation

We now study the system in the limit of large $\alpha_{\rm ann}$ values, without the assumption of KE. For this, we solve the non-integrated BE for the distribution function $f_{\chi}$, which is given by

$$\frac{\partial f_{\chi}}{\partial t}-Hp\frac{\partial f_{\chi}}{\partial p}=$$(23)

By change of variables from $\left(t,p_{\chi}\right)$ to $\left(x,q\right)$ such that $x=\frac{m_{\chi}}{T\left(t\right)}$ is the same as before and $q=a(x)p_{\chi}$ is the comoving momentum (with $a(x)$ the scale factor), Eq. (23) turns into a simple ordinary differential equation for each $q$,

$$xH\frac{d\bar{f}_{\chi,q}}{dx}=-c\left(x,q\right)\left(\bar{f}_{\chi,q}(x)-$$(24)

where

$$c\left(x,q\right)=\frac{2g_{\psi}m_{\psi}^{2}\alpha_{\rm decay}T}{g_{\chi}p_{$$(25)

with the notation $\bar{f}_{\chi,q}(x)=f_{\chi}\left(t(x),p_{\chi}\left(x,q\right)\right)$, $\bar{f}_{\chi,q}^{\rm eq}(x)\equiv\exp\left(-\frac{E_{\chi}\left(x,q\right)}{T$ and

$$\delta\equiv\frac{m_{\psi}^{2}-m_{\phi}^{2}-m_{\chi}^{2}}{2m_{\chi}^{2}}\,.$$(26)

The complete derivation can be found in Appendix B and in Refs. D’Agnolo et al. (2017); Garny et al. (2017).

Eq. (24) can be integrated in a similar way as the yield in the KE case in Eqs. (8) and (10):

$$\bar{f}_{\chi,q}(x)=e^{-\mathbf{A}_{\rm q}(x)}\left(\bar{f}_{q,\chi}\left(1$$(27)

with

$$\mathbf{a}_{\rm q}(x)\equiv\frac{c\left(x,q\right)}{xH},\quad\mathbf{b}_{\rm q$$(28)

and

$$\mathbf{A}_{\rm q}(x)\equiv\intop_{1}^{x}\mathbf{a}_{\rm q}\left(\xi\right)d$$(29)

Eqs. (28) and  (29) are analogous to Eqs. (9) and (11). We assume that at $x=1$ all the particles are at kinetic and chemical equilibrium with the SM bath.

Finally, the particle number density $n_{\chi}(x)$ is given by integration over all momenta, and the yield can be written as

$$Y_{\chi}(x)=\frac{1}{s\left(x=1\right)}\frac{g_{\chi}}{2\pi^{2}}\intop q^{2}$$(30)

The integral is dominated by the region of $q\sim m_{\chi}$. This can be understood by the following arguments:

• For small momenta ($q\ll m_{\chi}$) the distribution function at $x=1$ is similar ($\bar{f}_{\chi,q}^{\rm eq}\left(1\right)\sim\bar{f}_{\chi,m_{\chi}}^{\rm eq}$) and the collision rate is similar ($c\left(x,q\right)\sim c\left(x,m_{\chi}\right)$), but the contribution to the integral is small because of the $q^{2}dq$ measure.

• For high momenta ($q\gg m_{\chi}$) the initial value is exponentially smaller ($\bar{f}_{\chi,q}^{\rm eq}\left(1\right)\ll\bar{f}_{\chi,m_{\chi}}^{\rm eq}$), and the collision term exponentially increases with $q$.

Thus, the region $q\sim m_{\chi}$ is the most dominant one in the integral evaluation.

Figure 7: Impact of kinetic equilibrium. The ratio between the decay couplings obtained without and with kinetic equilibrium, as a function of DM mass. The various panels are for different mass splittings $\Delta$, while the curves inside each panel are obtained for different $\phi$ masses, represented by various values of $r=m_{\phi}/m_{\chi}$ relative to $r_{\rm max}=\Delta$ (see Eq. (5)). The calculations in this figure are made with constant $g_{*}$. Note that for $\Delta=0.5$, unitarity considerations typically limit the DM mass to $m_{\chi}\mathrel{\hbox to0.0pt{\lower 4.0pt\hbox{\hskip 1.0pt$\sim$}\hss}$ MeV, see Fig. 3.

IV.2 Numerical Evaluation of the Decay Coupling

Fig. 7 presents the ratio $\alpha_{\rm NKE}/\alpha_{\rm KE}$, namely the ratio between the couplings that reproduce the observed DM relic abundance without and with kinetic equilibrium, assuming equilibrium initial conditions, as a function of $m_{\chi}$, for various fixed values of $\Delta$ and $r$.

For all the curves with $\Delta\leq 0.2$, the NKE corrections to the coupling are relatively small (up to $20\%$). For $\Delta=0.5$, the corrections to the coupling might be up to an order of magnitude for high DM masses $m_{\chi}\sim 10^{7}-10^{8}$ MeV, but such high masses for this large mass splitting are unlikely in the mechanism due to the unitarity bound on the annihilation cross section (see Section III.1 and Fig. 3, which restricts to $m_{\chi}\mathrel{\hbox to0.0pt{\lower 4.0pt\hbox{\hskip 1.0pt$\sim$}\hss}$ MeV for large mass splittings). Moreover, the couplings ratio is close to unity around $1\,{\rm keV}$ and grows up to $10$ only after eleven orders of magnitude in DM mass for massless $\phi$ ($r=0$). For massive $\phi$ ($0<r<\Delta$), which is the case in the model we present latter in Section V, the slope in Fig. 7 decreases with the increase of $m_{\phi}$. We learn that the NKE corrections to the coupling values are rather small.

The following paragraphs are dedicated to explaining the proximity of the obtained couplings along with other features in Fig. 7. Since Eq. (27) has the same functional form as Eq. (10), comparison between the KE and NKE cases is discussed separately for the two INDY types of DM— the out of chemical equilibrium case which is relevant for small DM mass and small mass splittings; and the in chemical equilibrium case, relevant for large DM mass or large mass splittings.

IV.2.1 Small Mass Splitting Case

For low $\Delta$ values, the initial condition terms in Eqs. (12) and (27) dominate the relic abundance. To examine analytically such case and explaining why $\alpha_{{\rm KE}}\approx\alpha_{{\rm NKE}}$, we look at the ratio between the inverse decay rates $\mathbf{a}_{\rm q}$ of relevant $q$ mods ($q\sim m_{\chi}$) to the inverse decay rate of the yield in the KE case $\mathbf{a}$.

For the KE case, the inverse decay rate in the NR approximation is

$$xH\mathbf{a}(x)\approx\alpha_{\rm KE}\bar{\beta}m_{\psi}\frac{g_{\psi}}{g_{$$(31)

Assuming constant $g_{*}$ and $g_{*s}$ for simplicity (which provides $a=x$ and $q=p_{\chi}x$), the NKE inverse decay rate is

$$xH\mathbf{a}_{\rm q}(x)\approx\frac{2g_{\psi}m_{\psi}^{2}\alpha_{\rm NKE}}{g_{$$(32)

where we used NR approximations for the energy - $E_{\chi}\approx m_{\chi}+\frac{1}{2}\frac{p_{\chi}^{2}}{m_{\chi}}$ in the exponent, and $E_{\chi}\approx m_{\chi}$ in the denominator.

Since $\bar{\beta}\leq\bar{\beta}_{m_{\phi}=0}=\frac{m_{\psi}-m_{\chi}}{m_{\psi}}$, the term inside the hyperbolic sine function is small, and a linear approximation can be made, providing

$$xH\mathbf{a}_{\rm q}(x)\approx\frac{g_{\psi}}{g_{\chi}}\alpha_{\rm NKE}e^{-$$(34)

In the case of $m_{\phi}=0$, the ratio between the rates in Eqs. (34) and (31) is

$$\frac{\mathbf{a}_{\rm q}(x)}{\mathbf{a}(x)}\approx\frac{\alpha_{{\rm NKE}}}{$$(35)

and substituting $x_{*}=\Delta^{-1}$ provides

$$\frac{\mathbf{a}_{\rm q}\left(x_{*}\right)}{\mathbf{a}\left(x_{*}\right)}$$(36)

It follows from Eq. (36) that the ratio between the rates at $x\sim x_{*}$ is close to the ratio between the couplings $\frac{\mathbf{a}_{\rm q}\left(x_{*}\right)}{\mathbf{a}\left(x_{*}\right)}\sim$. Since this region of $x$ is contributes the most in evaluating the total probability to inverse decay, we have $\frac{\mathbf{A}_{\rm q}}{\mathbf{A}}\sim\frac{\alpha_{{\rm NKE}}}{\alpha_{{$.

Thus, the same probability for a particle to inverse decay is obtained when close values of decay coupling are taken, which provide the same (observed) relic abundance in both KE and NKE cases, as discussed above.

In the massive $\phi$ case, the ratio is the same as in Eq. (35), up to an additional factor of $\exp\left(\left(\frac{1}{2}x+\frac{1}{4x}\frac{q^{2}}{m_{\chi}^{2}}\right)$. For $x=x_{*}$, this factor equals

$$e^{\left(\frac{1}{2}x_{*}+\frac{1}{4x_{*}}\frac{q^{2}}{m_{\chi}^{2}}\right)$$(37)

Since $\frac{1}{\Delta}\frac{m_{\phi}}{m_{\chi}}<1$, we obtain $\frac{1}{\Delta}\left(\frac{m_{\phi}}{m_{\chi}}\right)^{2}<\Delta$ and $\frac{\Delta}{4}\frac{q^{2}}{m_{\chi}^{2}}\frac{m_{\phi}^{2}}{m_{\chi}^{2}}<$ and so the two exponents are small. Again, the ratio between the rates is close to the ratio between the couplings, and the same augment as in the massless $\phi$ case is valid, providing $\alpha_{\rm NKE}\sim\alpha_{\rm KE}$. Moreover, since in Eq. (37) the additional factor increases with $m_{\phi}$, $\alpha_{\rm NKE}$ should decrease to obtain the same rate, thus explaining why in Fig. 7 the $\Delta=0.05$ curves decrease slightly as $r$ increases.

As a final remark, we note that the maximum of $\mathbf{a}_{\rm q}$ is $q$ dependent and not exactly at $x_{*}$ , but the qualitative analysis presented above is not affected by the slight maximal rate time change.

IV.2.2 Large Mass Splitting Case

Here, we aim to understand why $\alpha_{{\rm NKE}}\sim\alpha_{{\rm KE}}$ even for relatively high $\Delta$ values. In this case, $\psi$’s decays dominate the relic abundance. Using the saddle point approximation, it can be shown that the phase space distribution function dependence on the mass scale (up to logarithmic corrections) is

$$\bar{f}_{\chi,q,\infty}\sim\left(\frac{\alpha_{\rm NKE}m_{\rm pl}}{q}\right)^{$$(38)

The complete derivation is given in Appendix A. Using Eq.  (30), the relic abundance’s dependence on $m_{\chi}$ is

$$Y_{\chi,\infty}\sim\left(\frac{\alpha_{\rm NKE}m_{\rm pl}}{m_{\chi}}\right)^{-$$(39)

and by matching the relic abundance in Eq. (39) with the observed value we produce the power-law relation

$$\alpha_{\rm NKE}\sim m_{\chi}^{1+\delta}.$$(40)

Finally, Comparing the NKE and KE cases (Eqs.  (19) and (40)) provides

$$\frac{\alpha_{\rm NKE}}{\alpha_{\rm KE}}\sim m_{\chi}^{\delta-\Delta}=m_{\chi}$$(41)

This remarkable result explains the low growth rate of the ratio for high $\Delta$ values observed in Fig. 7. For $\Delta=0.5,r=0$, the expected power is $\frac{1}{8}$, which is a small power by itself, and close to the slope in the figure, $\sim\frac{1}{11}$. For $\Delta=0.2$, the expected power reduces to an even smaller number $\frac{1}{50}$, explaining the growth of only $\sim 25\%$ after eleven orders of magnitude.

Additionally, increasing $r$ reduces the power, in correspondence with the $\Delta=0.5$ curves at Fig. 7. For $\Delta=0.1,0.2$, the power is so low that the leading term in Eq. (39) is not enough to explain the results, and a more subtle approach is needed.

In summary, we have verified that the results for $\alpha_{\rm decay}$ for INDY DM, which were obtained in the previous sections assuming kinetic equilibrium, are in close proximity to the obtained values in the absence of such an assumption, rendering our analysis valid.

V Model

The system considered in this work, described by the BEs (1), can be realized by a simple renormalizable toy model Frumkin et al. (2023a). We now present such a model and discuss its relevant parameter space which manifests the phases of INDY DM and coannihilation via decays, and discuss possible phenomenology.

We consider a dark sector with internal gauge symmetry $U(1)_{d}$ with gauge field $A_{d}$ (dark photon) and gauge coupling $e_{d}$ (dark charge). It contains two Dirac fermions $\chi$ and $\psi$, and a complex scalar $H_{d}$ with charges $0$, $1$, and $1$ respectively. Note that $\chi$ and $\psi$ take the same role as in the previous sections ($\psi$ decays to $\chi$), while $H_{d}$ does not explicitly appear in the mechanism as presented thus far. The general renormalizable Lagrangian contains mass terms for the fermions

$$\mathcal{L}_{\rm Dirac}=-m_{\chi}\bar{\chi}\chi-m_{\psi}\bar{\psi}\psi\,,$$(42)

spontaneous symmetry breaking (SSB) potential for the scalar

$$V\left(\phi\right)=-\frac{1}{2}\mu^{2}H_{d}^{*}H_{d}+\frac{1}{4}\lambda\left(H$$(43)

and a Yukawa interaction

$$\mathcal{L}_{\text{Yukawa}}=-yH_{d}^{*}\bar{\chi}\psi+h.c\,.$$(44)

The connection to the SM is through mixing between the dark photon and the photon ($A$) fields,

$$\mathcal{L}\supset-\frac{\epsilon}{2}F_{d}^{\mu\nu}F_{\mu\nu}\,.$$(45)

The field $H_{d}$ acquires a vacuum expectation value $v_{H_{d}}\equiv\left\langle H_{d}\right\rangle=\frac{\mu}{\sqrt{\lambda}}$ through SSB, providing a real scalar field $h_{d}$ with $\left\langle h_{d}\right\rangle=0$ and $m_{h_{d}}=\mu$, and mass term to the dark photon of Peskin and Schroeder (1995); Cheng and Li (1984)

$$m_{A_{d}}^{2}=2e_{d}^{2}v_{H_{d}}^{2}\,.$$(46)

Additionally, the Yukawa interaction term after SSB yields a mass between the fermions $\mathcal{L}\supset-yv_{H_{d}}\bar{\chi}\psi+h.c$. A rotation to the mass diagonalized fields $\chi^{\prime}$ and $\psi^{\prime}$, using the rotation angle $\vartheta$ defined by

$$\sin 2\vartheta=\frac{2yv_{H_{d}}}{\left(m_{\psi^{\prime}}-m_{\chi^{\prime}}$$(47)

removes the mass mixing. The term $\mathcal{L}\supset-e_{d}\bar{\psi}\not{A}_{d}\psi$ after the rotation yields an explicit interaction between $\psi^{\prime}$, $\chi^{\prime}$ and $A_{d}$:

$$\mathcal{L}\supset e_{d}\cos\vartheta\sin\vartheta\bar{\chi}^{\prime}\not{A}_{$$(48)

Since only the mass eigenstates ($\chi^{\prime}$ and $\psi^{\prime}$) are considered, in what follows we neglect the prime in denoting the fields.

V.1 Parameter restrictions

We impose several conditions on the parameters of the model, such that the decay $\psi\rightarrow\chi A_{d}$ and annihilation $2\psi\rightarrow 2A_{d}$ processes dominate the relic abundance.

First, for the decay to the dark photon to occur, the restriction $m_{\chi}+m_{A_{d}}<m_{\psi}$ must be imposed. Additionally, to prevent the decay to the dark Higgs $\psi\rightarrow\chi h_{d}$ from occurring (else it would dominate the relic abundance), we take $m_{\chi}+m_{h_{d}}>m_{\psi}$. The inequalities over the masses can be rewritten as a constraint over the dark charge

$$e_{d}<\sqrt{\frac{\lambda}{2}}\min\left\{\frac{m_{\psi}-m_{\chi}}{m_{h_{d}}}\;$$(49)

Next, for the dark sector to thermalize with the SM, the reaction rate $A_{d}\leftrightarrow e^{+}e^{-}$ must be fast enough along the evolution of the DM abundance. This provides a lower bound on $\epsilon$ which is imposed later in this section by requiring $\Gamma_{A_{d}\rightarrow e^{+}e^{-}}\geq H$ at $x=1$.

Finally, $h_{d}$ must follow equilibrium (to avoid co-scattering of $\chi$ with $h_{d}$ that is out of CE) and decay into lighter particles (otherwise it is a DM candidate by itself). This can be achieved in two ways. If $e_{d}<\frac{1}{2}\sqrt{\frac{\lambda}{2}}$ the reaction $h_{d}\rightarrow 2A_{d}$ is permitted and no further restriction is added. If $\frac{1}{2}\sqrt{\frac{\lambda}{2}}<e_{d}<\sqrt{\frac{\lambda}{2}}$, the reaction $h_{d}\rightarrow A_{d}e^{+}e^{-}$ is responsible for the dark Higgs decay. For $h_{d}$ to follow CE in the latter case, we demand $\Gamma_{h_{d}\rightarrow A_{d}e^{+}e^{-}}\geq H$ at $T=m_{h_{d}}$, which provides again a lower bound on $\epsilon$.

Figure 8: Model.Left: The dark charge $e_{d}$ verses the Yukawa coupling $y$, which reproduces the observed abundance for the model, for different DM masses. The lighter colored solid curves in the plot show the values obtained with decay and annihilation reactions only, without co-scattering. Here we set $\lambda=1$. Right: Allowed dark photon parameter space. Solid thick colored curves indicate the minimal kinetic mixing as a function of dark photon mass within the model, for various values of mass splitting $\Delta$ with $\lambda=1,\mu=2\left(m_{\psi}-m_{\chi}\right),e_{d}=0.32$. In this case, the constraint on $\epsilon$ comes from the thermalization of the sectors through the $A_{d}\leftrightarrow e^{+}e^{-}$ reaction. The dashed thick lines indicate the minimal kinetic mixing for different parameter choice: $\lambda=1,\mu=1.1\left(m_{\psi}-m_{\chi}\right),e_{d}=0.58$. Here, the constraint on $\epsilon$ comes from the $h_{d}\rightarrow A_{d}e^{+}e^{-}$ reaction rate (see Section V.1). In both cases the connection between the masses of the dark photon and the DM is $m_{A_{d}}=0.9\,\Delta m_{\chi}$. We show existing limits Fradette et al. (2014); Alexander et al. (2016); Chang et al. (2017); Hardy and Lasenby (2017); Pospelov and Tsai (2018); Banerjee et al. (2018); Aaij et al. (2018, 2020); Parker et al. (2018); Tsai et al. (2021) in shaded gray, with future projections Celentano (2014); Altmannshofer et al. (2019); Ilten et al. (2015); Alekhin et al. (2016); Ilten et al. (2016); Alexander et al. (2016); Caldwell et al. (2018); Berlin et al. (2018, 2019); Ariga et al. (2019); NA62 (2018); Tsai et al. (2021) indicated by the colored thin dashed curves.

V.2 Interactions

We now relate the model’s parameters to the couplings $\alpha_{\rm decay}$ and $\alpha_{\rm ann}$ in Eq. (2) and Eq. (3).

At the tree-level approximation, the decay coupling is given by Shtabovenko et al. (2020, 2016); Mertig et al. (1991)

$$\alpha_{\rm decay}=\frac{y^{2}\left(\left(m_{\chi}-m_{\psi}\right)^{2}-m_{A_{d$$(50)

Since the relevant values of $\alpha_{\rm decay}$ are very small (see Fig. 4), $y$ must be small as well.

By neglecting diagrams and terms with high orders of $y$, calculating the thermally averaged cross-section in the s-wave approximation provides Shtabovenko et al. (2020, 2016); Mertig et al. (1991)

$$\alpha_{\rm ann}^{2}=\frac{e_{d}^{4}\sqrt{1-\frac{m_{A_{d}}^{2}}{m_{\psi}^{2}}$$(51)

Additional reactions involving $\chi$ occur in the model, but their effect is small. Since $y$ is small, all the reactions with matrix element squared of order $\mathcal{O}\left(y^{3}\right)$ or above are insignificant, leaving only nine reaction types with $\overline{\left|\mathcal{M}\right|^{2}}\propto y^{2}$. In the following, we show that the affect of those nine reactions is small as well.

The fastest of the nine interactions is co-scattering with the dark photon $\chi A_{d}\rightarrow\psi A_{d}$, since it is proportional to the dark photon number density $n_{A_{d}}$. The ratio between its rate to inverse decay is

$$\frac{\text{Rate}_{\chi A_{d}\rightarrow\psi A_{d}}}{\text{Rate}_{\chi A_{d}$$(52)

In the parameter space considered in this work, the dark photon mass is settled near its threshold $m_{d}\lesssim m_{\psi}-m_{\chi}$ (see next subsection). From Eq.  (52), the ratio between the rates (and the co-scattering rate itself) decreases rapidly at $x\sim\frac{m_{\chi}}{m_{A_{d}}}$, which is near the temperature at which the inverse decay rate is maximal $x_{*}=\Delta^{-1}$ (see Eq. (14)). Thus, the inverse decays dominate the freeze-out process in such parameter choices, and co-scattering serves as a secondary process.

It is worth mentioning that since for $x<x_{*}$ the co-scattering rate might be higher, it can change the ‘effective’ initial conditions for the inverse decay process (a possibility mentioned in Garny et al. (2017)); however, we have not found significant change to the obtained coupling $y$ in the calculations done for this work. The small effect of initial conditions on $y$ is in correspondence with the results shown in Section III.2, which discusses the effect on changing the initial conditions on $\alpha_{\rm decay}$.

The remaining eight reactions have a lower rate compared to co-scattering with the dark photon. The rates of co-scattering involving the scalar $h_{d}$ ($\chi A_{d}\rightarrow\psi h_{d}$, $\chi h_{d}\rightarrow\psi h_{d}$ and $\chi h_{d}\rightarrow\psi A_{d}$) are slower than the rate of co-scattering with $A_{d}$ only by a factor of $\sim n_{h_{d}}^{\rm eq}/n_{A_{d}}^{\rm eq}$ or $\sim\left(n_{h_{d}}^{\rm eq}/n_{A_{d}}^{\rm eq}\right)^{2}$. The reactions with rate containing a factor of $n_{\psi}$ ($\chi\psi\rightarrow\psi\psi$ , $\chi\bar{\psi}\rightarrow\psi\bar{\psi}$, $\chi\bar{\psi}\rightarrow A_{d}A_{d}$, $\chi\bar{\psi}\rightarrow h_{d}A_{d}$ and $\chi\bar{\psi}\rightarrow h_{d}$) are exponentially suppressed as well ($\chi\bar{\psi}\rightarrow h_{d}$ is also forbidden in our parameter space choice).

V.3 Results

For the model to avoid cosmological constraints from $N_{\rm eff}$ and Big Bang Nucleosynthesis Parker et al. (2018); Banerjee et al. (2019); Lees et al. (2017); Andreas et al. (2012); Blümlein and Brunner (2014); Collaboration (2015); Tsai et al. (2022); Ibe et al. (2020); Sabti et al. (2020); Bauer et al. (2018), the parameters are chosen such that the dark photon (which is chosen to be the lightest particle) has mass of $m_{A_{d}}\geq 18\ {\rm MeV}$. High values of $\alpha_{\rm ann}$ are needed for high mass scales (see Fig. 3), which translates to a requirement of high $e_{d}$ values and large $m_{A_{d}}$ masses (see Eqs. (51) and (46)).

The left panel of Fig. 8 shows $y$ and $e_{d}$ values that reproduce the observed DM relic abundance for different DM masses, mass splittings and values of $r$. Co-scattering with the dark photon plays little role, as is evident from the proximity of the results with and without co-scattering. The analogy of Fig. 8 to Fig. 1 emphasizes how the model manifests the INDY and coannihilation via decay mechanisms, through Eqs. (51) and (50), for various parameters.

INDY DM can be expected to evade conventional searches such as direct and indirect detection of DM. The reason is that in the dark sector, all $\chi$ interactions are suppressed by factors of $y\sim 10^{-8}$, and for interaction with the SM, additional factors of $\epsilon\ll 1$ are present. For example, the obtained cross-section of $\chi e\rightarrow\chi e$ elastic scattering, which is given by Shtabovenko et al. (2020, 2016); Mertig et al. (1991); Peskin and Schroeder (1995)

$$\sigma_{e}=\frac{y^{4}\epsilon^{2}e^{2}m_{e}^{2}m_{\chi}^{2}}{4\pi e_{d}^{2}$$(53)

where $m_{e}$ and $e$ are the mass and electric charge of the electron, is smaller by several orders of magnitude than current reach of direct direction experiments Kuflik et al. (2017); Alexander et al. (2016); Essig et al. (2017); Liu et al. (2017); Essig et al. (2022). Note that while Eq. (53) is model dependent, the low values of $\alpha_{\rm decay}$ obtained in the general case (Fig. 4) suggests that $\chi$ interactions are expected to be weak in any model which realizes the mechanism.

The dark photon can, however, be searched for directly. The right panel of Fig. 8 shows the minimal value of $\epsilon$ needed for each $m_{A_{d}}$, for three choices of parameters within the model. Experimental constraints on a visibly decaying dark photon in the mass range $\mathcal{O}(10-1000)$ MeV currently constrain $\epsilon\lesssim{10^{-3}-10^{-4}}$ Fradette et al. (2014); Alexander et al. (2016); Chang et al. (2017); Hardy and Lasenby (2017); Pospelov and Tsai (2018); Banerjee et al. (2018); Aaij et al. (2018, 2020); Parker et al. (2018); Tsai et al. (2021), and most of the relevant parameter space is expected to be probed by future experiments, as demonstrated in the plot  Celentano (2014); Altmannshofer et al. (2019); Ilten et al. (2015); Alekhin et al. (2016); Ilten et al. (2016); Alexander et al. (2016); Caldwell et al. (2018); Berlin et al. (2018, 2019); Ariga et al. (2019); NA62 (2018); Tsai et al. (2021) (see Ref. Tsai et al. (2022) for a recent compilation of existing searches and future projections for visibly decaying dark photons). Additionally, the INDY dark matter candidate should be detectable through the direct production of $\psi$ particles that decay into $\chi$ and either additional invisible products or SM particles. We leave a detailed study of such a signal to future work.

VI Conclusions

In this work, we studied in detail the possibility that inverse decays play a meaningful role in setting the abundance of DM. We considered a framework of inverse decay interactions along with a self-annihilation process of the unstable particle, and examined the different phases which emerge in this setup. This includes a coannihilations via decay phase, in which the freezeout of the annihilation determines the relic abundance; an INDY phase, in which the freezeout of inverse decays determines the relic abundance; a FI phase, in which only the unstable particle decays account for the observed DM abundance; a FIFO phase, in which the DM abundance is given by freeze in and then freezeout of the DM; and a freezeout and decay phase, in which the decay plays a role only after the unstable particle freezes out.

We showed both numerically and by analytical approximations that when inverse decay reactions control the relic abundance of DM, much smaller couplings are required compared to other mechanisms such as the WIMP. (Interestingly, a chain of inverse decays can allow very heavy DM with larger couplings of ${\cal O}(1)$; for details see Ref. Frumkin et al. (2023b).) We further verified that relaxing the assumption of kinetic equilibrium does not impact (conceptually or numerically) the couplings or parameter space of INDY DM. Finally, we presented a renormalizable model that can realize the phases of DM described in this paper We mapped out the relevant parameter space within the model for the coannihilation via decays and INDY mechanisms presented in this work, and discussed the possible phenomenology for such DM. There is much room for further model-building efforts where inverse decays control the DM abundance. Interesting experimental signals at accelerators are anticipated from the long-lived particles of such theories, which could be used as a probe for inverse decaying dark sectors.

Appendix A Saddle Point Approximation for INDY DM In Chemical Equilibrium

In this Appendix, we show how we use the saddle point approximation to obtain Eqs. (18) and (38).

The relic abundance with kinetic equilibrium for an INDY DM that is in CE is dominated by the $\psi$’s decays term of Eq. (12),

$$Y_{\chi,\infty}=\intop_{1}^{\infty}e^{-\intop_{\eta}^{\infty}\mathbf{a}\left($$(54)

where the NR approximations for $\mathbf{a}(x)$ and $\mathbf{b}(x)$ can be written as

$$\mathbf{a}(x)=\mathbf{a}_{0}xe^{-\Delta x},\quad\mathbf{b}(x)=\mathbf{b}_{0}x^$$(55)

with

$$\mathbf{a}_{0} = \frac{3}{\pi}\sqrt{\frac{10}{g_{*}}}\frac{g_{\psi}}{g_{\chi}}\bar$$(56)

$$\mathbf{b}_{0} = \frac{3^{3}\cdot 5^{3/2}}{2^{2}\pi^{9/2}}\frac{g_{\psi}}{g_{*s}$$

Substituting the above into Eq. (54) yields

$$Y_{\chi,\infty}=\mathbf{b}_{0}\intop_{1}^{\infty}e^{-\frac{\mathbf{a}_{0}}{$$(57)

This integral can be approximated by the saddle point approximation,

$$\intop_{-\infty}^{\infty}e^{f(x)}d\eta\approx e^{f\left(\eta_{s}\right)}\sqrt{$$(58)

where we note $\eta_{s}$ as the distinct maxima of $f(x)$: $\frac{d}{d\eta}f\left(\eta_{s}\right)=0$. As long as $\eta_{s}>1$, which is the case for an INDY that is in CE, the lower bound of the integral in Eq. (57) can be extended from $1$ to $-\infty$ and Eq. (57) has the same form as Eq. (58) with

$$f(\eta)=-\frac{\mathbf{a}_{0}}{\Delta^{2}}e^{-\eta\Delta}\left(1+\eta\Delta$$(59)

The saddle point $\eta_{s}$ obeys

$$e^{-\left(1+\Delta\right)\eta_{s}}=\mathbf{a}_{0}^{-\frac{\left(1+\Delta\right$$(60)

and using Eqs. (57), (58), (59) and (60) one can approximate the relic abundance to

$$Y_{\chi,\infty}\approx\mathbf{b}_{0}\mathbf{a}_{0}^{-\frac{\left(1+\Delta$$(61)

Since the dependence of $\eta_{s}$ on $m_{\chi}$ is logarithmic (because $m_{\chi}$ appears only in $\mathbf{a}_{0}$ in Eq. (60)), the dependence of $Y_{\chi,\infty}$ on the DM mass (up to logarithmic corrections) is given by

$$Y_{\chi,\infty}\sim\mathbf{b}_{0}\mathbf{a}_{0}^{-\frac{\left(1+\Delta\right)}$$(62)

which is Eq. (18) of the main text.

The case without kinetic equilibrium is almost identical. Assuming the NR approximation $E_{\chi}\approx m_{\chi}$ (note that here we neglect the energy dependence on the momenta in the exponent of $\mathbf{a}_{\rm q}$, which is a more crude approximation than in Eq. (32)), one can write

$$\mathbf{a}_{\rm q}\left(x,q\right)\approx\mathbf{a}_{1}xe^{-\delta x},\quad$$(63)

with

$$\mathbf{a}_{1}=\frac{6}{\pi}\sqrt{\frac{10}{g_{*}}}\frac{g_{\psi}}{g_{\chi}}$$(64)

Substituting the definitions in Eq. (63) into Eq. (27) gives

$$\bar{f}_{\chi,q,\infty}=\mathbf{a}_{1}\intop_{x_{0}}^{\infty}e^{-\intop_{\eta}$$(65)

The saddle point $\eta_{s}$ in this case obeys

$$e^{-\left(\delta+1\right)\eta_{s}}=\mathbf{a}_{1}^{-\frac{\delta+1}{\delta}}$$(66)

and the approximations provide

$$\bar{f}_{\chi,q,\infty}=\mathbf{a}_{1}^{-\frac{1}{\delta}}\left(1+\delta-\frac$$(67)

Again, the $\eta_{s}$ dependence on the mass is logarithmic and the $\bar{f}_{\chi,q,\infty}$ dependence on the mass scale (up to logarithmic corrections) is given by the $\mathbf{a}_{1}$ term

$$\bar{f}_{\chi,q,\infty}\sim\left(\frac{\alpha_{\rm NKE}m_{\rm pl}}{q}\right)^{$$(68)

Using Eq. (68) and Eq. (30), it is possible to estimate the DM relic abundance in the NKE equilibruim case, as done in Eq. (39) in Section IV.

Appendix B Non-integrated Boltzmann Equation for $\chi$

In this Appendix, we derive Eq. (24) from Eq. (23). The approach used here is similar to that in Refs. D’Agnolo et al. (2017); Garny et al. (2017).

Simplifying the left hand side (LHS) of Eq. (23) is done by changes of variables, from $\left(t,p_{\chi}\right)$ to $\left(x,q\right)$, with a conversion defined by $t=t(x)$ and $p_{\chi}=q/a(x)$. Here, $q$ is the comoving momentum, and in these variables the LHS is

$${\rm LHS}=\left(\frac{\partial}{\partial t}-Hp_{\chi}\frac{\partial}{\partial p$$(69)

where we ignore the $\left(1+\frac{1}{3}\frac{d{\rm ln}g_{*s}(T)}{d\ln T}\right)$ factor afterwords in Section IV.

Simplifying the right hand side (RHS) of Eq. (23) is done in four steps. The first step is integrating over $\boldsymbol{p}_{\phi}$, using the delta function over the 3-momentum. We obtain

$${\rm RHS}=\frac{g_{\phi}}{2E_{\chi}}\intop\frac{1}{2E_{\phi}}d\Pi_{\psi}\left($$(70)

where $\left(\boldsymbol{p}_{\phi}\right)^{2}=\left(\boldsymbol{p}_{\psi}\right)^{2}+$.

The second step is integrating over $d\Omega_{\boldsymbol{p}_{\psi}}$, where the $z$-axis is chosen to be in the direction of $\boldsymbol{p}_{\chi}$. Since $\overline{\left|\mathcal{M}\right|^{2}}$ and the energies do not depend on $\phi_{\boldsymbol{p}_{\psi}}$, the integration over this angle is trivial and provides a factor of $2\pi$. Additionally, since $\overline{\left|\mathcal{M}\right|^{2}}$ does not depend on $\theta_{\boldsymbol{p}_{\psi}}$ the remaining delta function over the energy can be analytically integrated, obtaining

$${\rm RHS}=\frac{g_{\psi}g_{\phi}}{16\pi E_{\chi}}\intop_{p_{\psi}=0}^{p_{\psi}$$(71)

where $\Theta$ is the Heaviside function.

The third step is changing the remaining integration variable from $p_{\psi}$ to $E_{\psi}$. This provides

$${\rm RHS}=\frac{g_{\psi}g_{\phi}}{16\pi}\frac{1}{p_{\chi}E_{\chi}}\intop_{E_{$$(72)

where (it can be shown that $E_{\min}\geq m_{\psi}$)

$$E_{\min,\max}=\delta E_{\chi}\pm\frac{m_{\psi}^{2}\bar{\beta}p_{\chi}}{2m_{$$(73)

The fourth step is integrating over the energy $E_{\psi}$. Since $\overline{\left|\mathcal{M}\right|^{2}}_{\psi\rightarrow\phi\chi}$ does not depend on the energies of the particles and assuming both $\psi$ and $\phi$ are in chemical equilibrium in the Maxwell-Boltzmann distributions, Eq. (72) is simplified to

$${\rm RHS}=-c\left(x=\frac{m_{\chi}}{T},q=p_{\chi}a(x)\right)\left(f_{\chi}$$(74)

with $c\left(x,q\right)$ defined in Eq. (25) in the main text.

Note that there is no need to assume a Maxwell-Boltzmann distribution for both $\psi$ and $\phi$. Since the final result must be of the form $\frac{d\bar{f}_{q,\chi}}{dx}\propto-\left(\bar{f}_{q,\chi}-\bar{f}_{q,\chi}^{$, assuming the distribution of one particle determines the proper approximation for the second, such that the BE would be consistent.

Finally, the connection between the averaged matrix element squared and the decay rate is

$$\Gamma=\frac{1}{2m_{\psi}}g_{\chi}g_{\phi}\overline{\left|\mathcal{M}\right|^{$$(75)

and so in the language of the parameterization of Eq.  (2), we have

$$\overline{\left|\mathcal{M}\right|^{2}}_{\psi\rightarrow\chi\phi}=\frac{16\pi m$$(76)