New Pathways to the Relic Abundance of Vector-Portal Dark Matter

We fully explore the thermal freezeout histories of a vector-portal dark matter model, in the region of parameter space in which the ratio of masses of the dark photon $A^{\prime}$ and dark matter $\chi$ is in the range $1 \lesssim m_{A^{\prime}}/m_{\chi} \lesssim 2$. In this region $2 \rightarrow 2$ and $3 \rightarrow 2$ annihilation processes within the dark sector, as well as processes that transfer energy between the dark sector and the Standard Model, play important roles in controlling the thermal freezeout of the dark matter. We carefully track the temperatures of all species, relaxing the assumption of previous studies that the dark and Standard Model sectors remain in thermal equilibrium throughout dark matter freezeout. Our calculations reveal a rich set of novel pathways which lead to the observed relic density of dark matter, and we develop a simple analytic understanding of these different regimes. The viable parameter space in our model provides a target for future experiments searching for light (MeV-GeV) dark matter, and includes regions where the dark matter self-interaction cross section is large enough to affect the small-scale structure of galaxies.

We fully explore the thermal freezeout histories of a vector-portal dark matter model, in the region of parameter space in which the ratio of masses of the dark photon A and dark matter χ is in the range 1 m A /mχ 2. In this region 2 → 2 and 3 → 2 annihilation processes within the dark sector, as well as processes that transfer energy between the dark sector and the Standard Model, play important roles in controlling the thermal freezeout of the dark matter. We carefully track the temperatures of all species, relaxing the assumption of previous studies that the dark and Standard Model sectors remain in thermal equilibrium throughout dark matter freezeout. Our calculations reveal a rich set of novel pathways which lead to the observed relic density of dark matter, and we develop a simple analytic understanding of these different regimes. The viable parameter space in our model provides a target for future experiments searching for light (MeV-GeV) dark matter, and includes regions where the dark matter self-interaction cross section is large enough to affect the small-scale structure of galaxies. The particle nature of dark matter (DM) remains a mystery, whose solution requires us to search beyond the Standard Model (SM). There are a great many suggestions for new physics particles that might solve the DM puzzle. One well-studied class of DM candidates is the weakly-interacting massive particles (WIMPs). WIMPs are theoretically attractive because they naturally arise in various Beyond-Standard Model (BSM) theories of new weak-scale physics, and because the thermal production of WIMPs, through their 2 → 2 annihilations to SM particles, naturally leads to the correct DM relic abundance. However, with increasingly strong experimental constraints being placed on the WIMP scenario, we are also motivated to consider alternative scenarios where other interactions control the final DM abundance.

CONTENTS
There has been considerable recent interest in exploring thermal relic scenarios that naturally produce DM at light (sub-GeV) masses (see for example Ref. [1]), as existing direct detection constraints are much less sensitive to sub-GeV mass DM (e.g. Refs. [2][3][4][5]). Existing beam dump experiments are sensitive to sub-GeV DM but leave much of the parameter space unconstrained [6][7][8][9][10][11][12][13][14][15]. New accelerator and direct detection experiments will soon explore the parameter space of light DM with unprecedented sensitivity (see Ref. [1] and references therein); consequently, it is important to understand the landscape of models which naturally populate this sub-GeV region.
Previous studies have identified a mechanism for thermally producing sub-GeV DM in which strong 3 → 2 self-annihilations among DM particles control the thermal relic abundance. This strongly-interacting-massiveparticle (SIMP) scenario naturally leads to stronglycoupled DM (α D ∼ 1) with mass similar to the QCD scale (m χ ∼ 10-100 MeV) [16]. The natural emergence of the strong scale in the thermal SIMP scenario makes it a particularly attractive framework. In this scenario, the DM and SM sectors remain in thermal equilibrium throughout freezeout via elastic scattering between DM and SM particles.
An alternative thermal production mechanism for light DM arises when this condition is relaxed; in the elastically decoupling relic (ELDER) scenario the DM and SM sectors thermally decouple through the elastic DM-SM scattering while strong 3 → 2 self-annihilations are still active [17,18]. In the ELDER scenario, although thermal freezeout proceeds through the 3 → 2 DM selfannihilations, the DM relic abundance is nevertheless determined by the decoupling of DM-SM elastic scattering. This is achieved through a dark sector process called "cannibalization" [19], which occurs immediately after elastic decoupling and proceeds until 3 → 2 freezeout. During cannibalization, while the DM and SM sectors are thermally secluded, 3 → 2 DM self-annihilations convert mass to kinetic energy and heat the dark sector. As a result, the dark sector temperature evolves slowly (logarithmically as a function of SM temperature) during cannibalization, and likewise, the DM abundance evolves slowly. This leads to a DM relic abundance that is primarily determined by its value at kinetic decoupling. The ELDER scenario also naturally leads to MeV-GeV mass DM.
Distinctive thermal production mechanisms for light DM have also been realized in the well-studied vectorportal DM model of a Dirac fermion DM particle χ charged under a hidden U(1) gauge symmetry with dark gauge boson A , which is coupled to the SM photon through kinetic mixing. In the region of parameter space in which the dark photon is more massive than the DM (r ≡ m A /m χ > 1), the kinematically suppressed 2 → 2 annihilations of DM to heavier A s (χχ → A A ) can control the relic abundance. In this "forbidden DM" (FDM) mechanism [20,21] the exponential suppression of the 2 → 2 process setting the relic abundance of DM naturally gives rise to DM exponentially lighter than the weak scale. The FDM mechanism was shown to be a viable mechanism for producing sub-GeV DM.
More recently, Ref. [22] showed that in the region of parameter space of the dark photon model in which 1.5 r 2, the kinematic suppression of the χχ → A A annihilation process is compensated for by a kinematically allowed 3 → 2 (χχχ → χA ) annihilation channel, which can then play a dominant role in setting the thermal relic abundance of DM. This "not-forbidden dark matter" (NFDM) scenario is analogous to the thermal SIMP scenario in that 3 → 2 processes can determine the DM relic abundance, realized in the simple and wellstudied vector-portal DM model. The NFDM scenario was also demonstrated to be a viable mechanism for naturally producing sub-GeV DM. In both the FDM and NFDM scenarios, the DM and SM sectors were assumed to remain thermally coupled throughout the freezeout of DM.
In this paper, we extend both of these frameworks to consider cases in which the DM and SM sectors are allowed to kinetically decouple during thermal freezeout of the DM. We fully explore the 1 r 2 region of parameter space of the dark photon model, in which the kinematically suppressed 2 → 2 (χχ → A A ) channel and the kinematically allowed 3 → 2 (χχχ → χA ) channel play important roles in controlling thermal freezeout, and relax the condition that kinetic equilibrium is maintained between the two sectors throughout the freezeout process. We find a rich set of novel cosmological histories leading to a range of different mechanisms for obtaining the correct DM relic density. Among these, we identify a general class of mechanisms in which the DM relic abundance is determined by processes controlling the kinetic decoupling of the DM and SM sectors (which we call the KINetically DEcoupling Relic -KINDER). This KINDER scenario in the dark photon model generalizes the ELDER scenario to cases in which multiple processes control the thermal coupling between dark and SM sectors, and in which a 3 → 2 annihilation process among multiple dark sector species supports heating of the dark A 0 sector.
The outline of our paper is as follows. In Section II, we describe the dark photon model in the 1 r 2 region we consider, including the primary interactions controlling chemical equilibrium in the dark sector, and those between the dark sector and SM particles. In Section III we discuss general features of dark sector freezeout in our model, including the relevant interaction processes, the Boltzmann equations, which describe the thermodynamic evolution of the system, and the freezeout conditions of relevant processes. In this section, we also classify three thermodynamic phases (A, B, and C), which generally describe the various stages of the thermal histories realized in our model.
In Sections IV and V, we characterize the thermal freezeout histories for 1.5 r 2 and 1 r 1.5 respectively. In each case, we identify a rich set of freezeout histories and analytically determine the parameter space regions where they occur. These different histories are naturally classified into specific regions in the -α D plane, where describes the mixing between the dark photon and the SM photon, and α D is the dark sector coupling. In Section IV we study the 1.5 r 2 region of our model, where the 2 ↔ 2 process freezes out before the 3 ↔ 2 process; the possible histories can be classified into the WIMP, NFDM and KINDER regimes. In Section V, we examine the 1 r 1.5 region of our model, where the 3 → 2 process freezes out prior to the 2 → 2 process, and find four distinct regimes in addition to the WIMP regime (Regimes I -IV). In Section VI we discuss the relevant experimental and cosmological constraints; finally, in Section VII we summarize our conclusions.
Throughout this paper, we make use of Planck 2018 cosmological parameters [23], using the TT,TE,EE+lowE+lensing results; we take the DM abundance to be the central value of Ω χ h 2 = 0.12, with h = 0.6736. All quantities are expressed in natural units, with = c = k B = 1. Finally, we use many different symbols for approximations in this paper, and have attempted to keep them consistent with the following definitions: (i) we use " " when the approximation is a physical limit, e.g., a nonrelativistic limit; (ii) we use "≈" for statements that are true within an order of magnitude, but which we will take to be an equality for the purpose of analytic results; (iii) finally, we use "∼" for statements that are true within an order of magnitude, but we do not use the fact either analytically or numerically.

II. MODEL
In the mass basis, the Lagrangian of the dark photon model we consider is: where the gauge coupling is α D = g 2 D /4π, and / D ≡ / ∂ −ig D / A . The dark photon A kinetically mixes with the SM photon, giving rise to a small coupling between the dark photon and the SM electromagnetic current J µ EM , set by the kinetic mixing parameter . The value of can naturally range from as small as 10 −13 up to 10 −1 [24]. The hidden U(1) symmetry can be spontaneously broken through a Higgs-like mechanism with the dark Higgs taken to be heavy enough to be excluded from this lowenergy effective description, since we will always be considering energies m χ , m A . The kinetic mixing generates the tree-level interactions between dark and SM particles shown in Fig. 1 We are primarily interested in scenarios in which the dominant DM-number-changing interactions are the χχχ ↔ χA (3 ↔ 2) process and the kinematically suppressed χχ ↔ A A (2 ↔ 2) process, shown in Fig. 2. This restricts us to the region of parameter space in which 1 r 2.
At lower values of r, where r < 1, the dominant process controlling thermal freezeout is χχ → A A (which is then kinematically allowed), and the A decays promptly to SM particles. This regime is strongly ruled out by cosmic microwave background (CMB) constraints on the annihilation cross section of DM into SM particles for m χ 10 GeV [25].
At higher values of r, where r > 2, the s-channel annihilation of χχ → ff via an off-shell A dominates the DM-number-changing interactions: the χχ → A A process is very kinematically suppressed and the χχχ → A χ process reduces to a scattering process among the χs as the final-state A promptly decays back to χχ. Dark sector freezeout proceeds via the classic WIMP freezeout scenario, which also runs into stringent CMB constraints on the s-wave annihilation of Dirac fermion dark matter below 10 GeV [23].
In the intermediate (1 r 2) region of interest to us, which of the 2 → 2 or 3 → 2 processes dominates during thermal freeze-out depends on the ratio r. The 2 → 2 process receives a kinematic suppression from χ particles annihilating into heavier A particles (with an exponential factor of the form e −2(r−1)mχ/T ), while the 3 → 2 receives a Boltzmann suppression from an extra factor of χ number density in the initial state (with an exponential factor of the form e −mχ/T ). In the lower half of the range in r we consider (1 < r 1.5), the Boltzmann suppression is more severe for the 3 → 2 process, and therefore the 2 → 2 process dominates during thermal freezeout. In this regime, and for the case in which the DM and SM sectors remain thermally coupled throughout thermal freezeout, the 2 → 2 process determines the relic abundance; this is the FDM scenario described in the introduction.
In the upper half of the range in r we consider (1.5 r < 2), in contrast, the large kinematic suppression of the 2 → 2 process renders it subdominant to the 3 → 2 process at freezeout. In this regime, and for the case in which the DM and SM sectors remain thermally coupled throughout thermal freezeout, the 3 → 2 process determines the thermal relic abundance; this is the NFDM scenario described in the introduction.

III. DARK SECTOR FREEZEOUT
Before we detail all of the different regimes in which the dark sector can evolve to obtain the final DM relic density, we will begin by discussing some general features of the dark sector freezeout in this model. By "dark sector freezeout", we mean the cosmological evolution from the initial state, when both sectors are in thermal equilibrium, to the point where the DM has attained its final comoving relic abundance.

A. Thermodynamic Variables
Throughout freezeout, for the parameter space we consider, χ and A remain in thermal equilibrium with each other through χ-A scattering. The dark sector can therefore be described by a single dark sector temperature T . The general expressions for the number densities of the particles in the nonrelativistic limit are: where g χ = 2 and g A = 3 are the numbers of degrees of freedom associated with each particle. The factor of two in Eq. (2) accounts for the fact that we are including both χ and χ in the definition of n χ . We have also included effective chemical potentials µ χ and µ A which are in general nonzero; we denote the number densities of χ and A with zero chemical potential as n χ,0 (T ) and n A ,0 (T ) respectively. We will also frequently use the inverse dimensionless temperatures, x ≡ m χ /T and x ≡ m χ /T . The energy densities and pressures of χ and A are related to their number densities by: and where the Maxwell-Boltzmann distributions with zero chemical potential for ρ χ,0 (T ) and P χ,0 (T ) are given by Similar relations hold for A . The entropy of the dark sector is conserved when no heat is transferred between the dark sector and the SM through processes that involve both dark sector and SM particles. The entropy density of the dark sector is Entropy conservation of the dark sector in the limit where heat transfer processes are inefficient is a useful fact that we will use extensively in obtaining an analytic understanding of our results. When the dark sector entropy is conserved, d(s D a 3 )/dt = 0, where a is the expansion scale factor. Since we will be discussing the time evolution of n χ and n A frequently in the context of analytic estimates, we will derive here several expressions related toṅ χ andṅ A that will be useful throughout the paper. First, taking the time derivative of n χ giveṡ where we have used dT /dt −HT . We will often make the approximation that m χ /T 1 during freezeout, and so the term 3/(2T ) can often be neglected, unless µ χ ∼ m χ . Similarly, We will often be interested in comparing the final number density of the dark matter after the dark sector completely freezes out to the number density required to achieve the relic abundance of dark matter today. Defining Y χ ≡ n χ /s SM , where s SM is the entropy density of the SM sector after the dark sector has completely decoupled, the correct relic abundance is obtained when [26] Y χ ≡ n χ s SM = 4.32 × 10 −10 GeV m χ .

B. Relevant Processes
In the conventional WIMP regime, DM freezes out through the process χχ → f f , where f is a SM fermion. Once the mixing parameter 10 −5 -10 −4 , however, χχ → f f freezes out while other dark sector processes are still active, and these processes play a significant role in the freezeout of the dark sector [27,28].
Outside the WIMP regime, there are four main processes that play important roles during the freezeout of the dark sector when 1 r 2: 1. The 2 ↔ 2 dark sector process, χχ ↔ A A . This process was shown in Ref. [27] to be responsible for the freezeout of the dark sector for 1 r 1.5, under the assumption that the dark sector was in full thermal equilibrium with the SM. As we described r  I. List of f (r) and g(r) values, as defined in Eqs. (12) and (13), evaluated at typical r-values of interest in this paper.
in the introduction, this process is kinematically forbidden for r > 1 for stationary χ particles, leading to a velocity-averaged annihilation cross section that is exponentially suppressed as a function of the dark sector temperature T . Explicitly, the annihilation cross section is given by [27] σv χχ→A A = n 2 We provide the expression for σv A A →χχ in App. A; to make our analytic estimates more convenient, however, we parametrize this annihilation cross section as follows: where g(r) is a function of r that captures the nontrivial r-dependence. Typical values of g(r) are shown in Table I.
As r increases, the rate of the forward process becomes exponentially more suppressed as the mass difference between χ and A increases. Note that in the forward direction, χχ → A A removes kinetic energy from the dark sector; the rate of the forward reaction also becomes exponentially suppressed as T decreases, since less kinetic energy is available to χ particles for conversion into the rest mass of A particles.
2. The 3 ↔ 2 dark sector process, χχχ ↔ A χ. For 1.5 r 2, the freezeout of the dark sector is mainly controlled by this process, as examined in Ref. [28], once again under the assumption of a dark sector in thermal equilibrium with the SM. The forward process is a 3 → 2 process, with velocityaveraged annihilation cross section given by where f (r) encodes the nontrivial r-dependence of the cross section; once again, the full expression for σv 2 χχχ→χA is given in App. A. For ease of notation, we will drop the subscript on the thermally averaged cross section from here on, unless it is needed to avoid ambiguity. Typical values for f (r) across the range of r considered in this paper are shown in Table I. Note that the forward reaction converts rest mass to kinetic energy, and heats the dark sector, similar to other 3 → 2 processes found in cannibal dark matter models [17,18,[29][30][31]. 3. A ↔ f f . The dark photon kinetically mixes with the SM photon, and can decay into a pair of SM fermions. This process is an important numberchanging process for A particles, and is one of two important processes responsible for transferring energy between the two sectors. The decay width Γ of A is given in full in App. A.
4. χf ↔ χf . This elastic scattering process, and all possible processes related by conjugation, allows χ to directly transfer energy to or from the SM. This process as well as A ↔ f f together determine how efficiently energy gets transferred between the two sectors. Once both χf ↔ χf and A ↔ f f become sufficiently inefficient, the dark sector and the SM can lose thermal contact and kinetically decouple, falling out of thermal equilibrium.
There are additional 3 ↔ 2 dark-sector-only processes that we do not consider, such as χχA → A A and A A A → χχ. Since we are only considering m A > m χ , these 3 → 2 processes have rates that are parametrically suppressed by at least one power of exp(−(r − 1)m χ /T ) compared to χχχ → A χ, and α D times at least one power of exp(−m A /T ) compared to χχ → A A . These slower processes are therefore relatively unimportant compared to the much faster 3 ↔ 2 and 2 ↔ 2 processes shown here.
We also neglect the processes A f ↔ γf and A γ ↔ f f : these processes are suppressed by an additional factor of the electromagnetic fine structure constant α EM relative to A → f f , and are also Boltzmann suppressed by n A n χ relative to χf → χf . Consequently, they never control when thermal decoupling between the two sectors occurs. They also do not play any important role based on the analytic understanding that we will develop below; they may only appear as terms proportional to n A − n A ,0 (T ) in the Boltzmann equations, and can therefore be treated as small corrections to energy transfer rate arising from decays.

C. Boltzmann Equations
The evolution of the system is governed by the coupled Boltzmann equations for the number densities of χ and A , n χ and n A , respectively, along with their energy densities ρ χ , ρ A and pressures P χ , P A : and where n f is the number density of charged SM particles, which for simplicity we assume to consist only of electrons and positrons. This assumption is justified because we are considering sub-GeV dark matter and dark photons, so thermal equilibrium between the SM and dark sector typically holds down to temperatures of T 100 MeV, at which point all other SM particles have annihilated or decayed away. Note that all dark sector (SM) variables are evaluated at the dark sector temperature T (SM temperature T ) unless otherwise stated. The prefactors for each term account for our convention of including both χ and χ in n χ , and for initial state symmetry factors. Our convention, as well as the derivation of the dark sector annihilation cross sections, can be found in Ref. [28]. We take the limit of nonrelativistic χ and A for the A → ff and χχ → ff energy transfer rates. Details on the energy transfer rate for elastic scattering χf → χf can be found in App. B; in particular, we highlight the fact that we have calculated σvδE χf →χf analytically without assuming that f is relativistic, which is to our knowledge a new result. This result is important when Eqs. (14)-(16) contain three unknowns: n χ , n A and T , and can be solved numerically for the coupled evolution of these variables as a function of the SM temperature T . The numerical solution of these equations is used for all of the results throughout the paper.
We will also rely significantly on analytic approximations to gain some intuition for these results. To this end, it is useful to write the energy density Boltzmann equation Eq. (16) in the nonrelativistic limit. Expanding the energy densities to first order in 1/x , which is a small parameter once T, T m χ , we find ρ χ m χ n χ 1 + 3 2x , and similarly for ρ A . We have also made the approximation dx/dt Hx. With these expansions, we obtain We have neglected DM annihilation into SM fermions in this analytic estimate for simplicity, since this process is typically not important in the regions of parameter space we will be interested in. We also find numerically that (x/2x 2 )dx /dx O(1) in all scenarios, and thus can be neglected in Eq. (18). The simplified Boltzmann energy density equation to leading order then readṡ Comparing this with the sum of the χ and A number density Boltzmann equations, Eqs. (14) and (15), which is given bẏ we finally obtain the following compact expression for the χ number density evolution: Comparing this expression with the number density Boltzmann equation for χ, we find With this relation, we can also reformulate the number density Boltzmann equation for A aṡ These equations show that in the nonrelativistic limit, the Boltzmann equations establish certain relations between the rates of the various processes, determined ultimately by number and energy conservation. These equations will prove to be extremely useful for gaining analytic understanding of our numerical results.

D. Fast Reactions and Freezeout
To gain an understanding of the freezeout behavior of our dark sector, it is useful to understand when processes are occurring at rates fast enough to influence the freezeout process, and when they cease to be important. For temperatures T m χ , the rates of all of the process are generally fast, i.e. the rates of all processes in one direction are all much larger than the Hubble rate. For example, the 3 ↔ 2 process is considered fast when While a process is fast, the corresponding terms in square brackets in Eqs. (14) and (15) will generically be small, e.g. for the 3 ↔ 2 process, such that otherwise, the 3 ↔ 2 process can change the number densities of both χ and A within a time much faster than the Hubble time, until Eq. (26) is satisfied. Similarly, the 2 ↔ 2 process is fast when with such that Once T m χ , the number densities of both χ and A are Boltzmann suppressed and rapidly decrease. At some point, the forward rates of these processes become comparable to the Hubble rate, and the process freezes out. For the 3 → 2 process, this happens when and for the 2 ↔ 2 process, Similar results hold for χχ → f f , just like in the conventional WIMP scenario. The approximate relations found in Eqs. (25) and (28) when the 3 ↔ 2 and 2 ↔ 2 processes are fast can be rewritten in terms of the effective chemical potential µ χ and µ A as and respectively. Note that when both processes are fast, these relations together enforce µ χ ≈ µ A ≈ 0. For processes that are responsible for transferring heat between the SM and dark sector, the criterion for when these processes are "fast" depend on how much heat is generated/removed due to the 2 ↔ 2 and 3 ↔ 2 processes described above. Since the energy density of the dark sector for T m χ is dominated by the χ as n χ n A , the rate of change of dark sector energy density per dark sector particle is given approximately by m χṅχ /n χ ; processes are considered "fast" if they can transfer heat between the sectors at a comparable rate.
As discussed in Sec. III B, the two most important processes transferring energy between the two sectors are A ↔ f f and χf ↔ χf . Let us first focus on the process A ↔ f f . In scenarios where both the 3 ↔ 2 and 2 ↔ 2 processes are fast, the number densities of the dark sector particles are given by n χ,0 (T ) and n A ,0 (T ). When T m χ , A ↔ f f is generally fast enough to maintain thermal equilibrium between the two sectors, so that T = T . However, once m χ > T , n A ,0 (T ) drops rapidly, and the number densities of the dark sector n χ,0 (T ) and n A ,0 (T ) evolve to a point where After this point, the term on the left-hand side starts to become small relative to the right-hand side, and A ↔ f f becomes ineffective at maintaining both sectors in thermal equilibrium. Similarly, the dark sector number densities can evolve to a point where after which χf ↔ χf is too slow to maintain thermal equilibrium. Once both Eq. (34) and (35) have been met, kinetic decoupling occurs, and the dark sector temperature T starts to diverge from the SM temperature T . Keep in mind that σvδE χf →χf is proportional to (T − T )/T ; the comparison made in Eq. (35) is therefore between the heat transfer rate when |T − T |/T ∼ O(1) and the energy lost due to n χ,0 decreasing. We are now ready to understand the broad features of the thermodynamic evolution of the dark sector. There are three thermodynamic phases that the dark sector in our model may go through: 1. Thermodynamic phase A: dark sector in thermal equilibrium with the SM. Interactions between the dark sector and the SM allow the two sectors to exchange heat. If these interactions are sufficiently fast, the dark sector stays in thermal equilibrium with the SM with T = T , and the number densities of χ and A are simply given by n χ,0 (T ) and n A ,0 (T ); 2. Thermodynamic phase B: T = T with zero chemical potential.
Once A → f f and χf → χf become too slow, the dark sector kinetically decouples, and develops a temperature different from T . The 2 ↔ 2 and 3 ↔ 2 dark sector processes can inject or remove heat from the dark sector. While both processes are fast, Eqs. (32) and (33) enforce 3. Thermodynamic phase C: T = T , with nonzero chemical potential. If either the 3 → 2 or the 2 → 2 process freezes out after the SM-dark sector processes become slow, χ and A develop a chemical potential µ χ (T ) and µ A (T ) respectively, according to either Eqs. (32) and (33).
In some parts of parameter space in the models we study, the dark sector goes through all three phases sequentially; in other parts of parameter space, a nonzero chemical potential develops once T starts diverging from T , leading to a direct transition from phase A to C without spending any significant time in phase B.
Previous studies investigating this model [27,28] have assumed that the dark sector only stays in thermodynamic phase A, with Ref. [27] making the further assumption that n A = n A ,0 (T ) throughout in their thermally coupled model. However, we shall see that for values of as large as 10 −5 , the dark sector does not stay in thermodynamic phase A throughout the process of freezeout, changing the dependence of the relic abundance on the model parameters drastically.
Throughout this paper, we will mostly be interested in values that are small, of order 10 −5 or smaller. However, if is too small, the dark sector and the SM sector need not have been in thermal contact at any point, calling into question the basic assumption we make that the two sectors start out in thermal equilibrium. To obtain an estimate for the minimum value of above which we are guaranteed to have the dark sectors in thermal equilibrium at T ∼ m χ , we follow Ref. [32], and set this minimum value of to be when the f f → A rate exceeds the Hubble rate at T = m A . When this condition is met, A particles can be produced at a rate much faster than Hubble at T ∼ m χ , allowing the whole dark sector to come into chemical equilibrium with the SM prior to the onset of the Boltzmann suppression from A and χ going nonrelativistic. This condition can be written as [32] where ζ is the Riemann zeta function, and M pl is the Planck mass. Using the expression for Γ in App. A and setting T = m A , we obtain the following estimate for eq , the minimum value of at which thermal equilibrium is guaranteed by T ∼ m χ : In practice, experimental constraints will limit us to values of 10 −8 ; we can therefore safely assume the dark sector to be thermally coupled to the SM at T ∼ m χ throughout this paper.

IV. 1.5 r 2
We begin our discussion of the freezeout of the vectorportal dark matter model with 1.5 r 2. For these values of r, the 2 ↔ 2 process freezes out before the 3 ↔ 2 process. Under the assumption that the dark sector stays in thermodynamic phase A with T = T , this regime -which we call the "classic not-forbidden dark matter (NFDM)" regime -was studied in Ref. [28], and was found to be a viable model for sub-GeV dark matter with appreciable self-interaction rates and thus the potential to affect the small-scale structure of galaxies. Here, we explore 1.5 r 2 including the temperature evolution of the dark sector.

A. "Classic Not-Forbidden" Regime
For sufficiently small values of with r 1.5, the 3 ↔ 2 process eventually freezes out later than χχ → f f -the process that controls conventional WIMP freezeoutand starts to become the main process that controls the final abundance of χ. This transition occurs when i.e. when both processes freeze out at roughly the same time. Using the analytic expressions for the quantities above, we obtain an estimate for N/W , the value of that sets the boundary between the 'classic NFDM' regime and the WIMP regime: where x f ≡ m χ /T f , and T f is the temperature at which freezeout of either of these two processes occur. g * is the effective number of relativistic degrees of freedom that enters into the Hubble parameter, H(T ) = 1.66g Further requiring that the final relic abundance of DM is equal to the observed one today gives a relation between α D , and m χ . In the WIMP regime, where freezeout is controlled by χχ → f f , the correct relic abundance is obtained when Eq. (10) is satisfied. This allows us to predict: as the boundary between the conventional WIMP-like regime and the "classic NFDM" regime when the correct relic abundance is achieved. For < N/W , the freezeout of the 3 ↔ 2 process determines the abundance of DM, and the parameters that generate the correct relic abundance become virtually independent of , provided that is large enough that the system remains in thermodynamic phase A (i.e. T = T ) throughout freezeout.

B. Kinetically Decoupling Relic (KINDER) Regime
As decreases further, processes that exchange energy between the dark and SM sectors become gradually less efficient; eventually, thermal equilibrium between the two sectors is lost even prior to 2 ↔ 2 freezeout. This scenario, which we call the kinetic decoupling relic (KINDER) regime, is starkly different from the "classic NFDM" regime explained above. Notably, the abundance of DM after freezeout is governed primarily by when kinetic decoupling occurs, and therefore depends on both and α D . With thermal equilibrium between the two sectors lost prior to the freezeout of dark sector processes, the dark sector now goes through the different thermodynamic phases described in Sec. III D.

General Features
In Fig. 3, we show the abundances of χ and A , as well as the dark sector temperature T as a function of x for our benchmark parameter values in the KINDER regime: m χ = 10 MeV, α D = 1, = 4 × 10 −8 , and r = 1.8. For ease of presentation, we plot the abundance as m χ Y χ and m χ Y A , where Y i is defined in Eq. (10). In Fig. 4, we show the number density and energy density rates for the relevant dark sector processes; explicitly, these are the terms for each process that appear on the right-hand side of Eqs. (14) and (15) divided by n χ for number density rates, and the right-hand side of Eq. (16) divided by n χ for energy density rates. At this parameter point (which is representative of the KINDER regime), the dark sector freezeout proceeds through the following stages: 1. Kinetic decoupling, transition from thermodynamic phase A to B. While either χf → χf and A → f f occur at rates larger than or comparable to the kinetic energy production rate of χ (i.e. the left-hand side of Eqs. (34) or (35) are large compared to the RHS, kinetic equilibrium between the dark sector and SM particles is maintained at a common temperature T = T . Once this is no longer true, i.e. after both χf → χf and A → f f become slow, kinetic decoupling occurs, and T begins to diverge from T . For our benchmark parameter values, kinetic decoupling occurs when A → f f becomes slow, as shown in Fig. 4. tinues until the 3 ↔ 2 process freezes out. This is an extension of the conventional cannibal dark matter scenario that we will investigate in greater detail below. 4. Freezeout of 3 ↔ 2 process. Finally, the 3 → 2 rate falls below the Hubble rate. With no other active number changing processes, the dark matter number density n χ evolves proportionally to a −3 .
Because the slow evolution of Y χ takes place from the time of kinetic decoupling until the freezeout of the 3 ↔ 2 process, the DM thermal relic density is governed mainly by the kinetic decoupling process. In this regime, the vector-portal DM model therefore shares many similarities with elastically decoupling (ELDER) dark matter [17], with the main differences being the existence of thermodynamic phase C mentioned in the last paragraph, and the fact that kinetic decoupling in vector-portal DM is frequently governed by A ↔ f f , instead of elastic scattering processes, i.e. χf ↔ χf . The dark sector entropy is also not fully conserved due to the existence of Similarly to the boundary between the WIMP and "classic NFDM" regimes, we can estimate the value of at which we transition from the KINDER regime to the "classic NFDM" regime, by finding the value of for which kinetic decoupling and 3 ↔ 2 freezeout occur at roughly the same time. We find that A ↔ f f is often the process that governs kinetic decoupling, and so the boundary between these regimes occurs at the value of = K/N where both Eqs. (30) and (34) are satisfied at the same SM temperature T . Analytically, we find K/N ∼ 10 −7 e 9.9(r−1.6) α D 1.0 as the boundary between the 'classic NFDM' regime and the KINDER regime, with x f denoting the dimensionless inverse temperature at the freezeout of the 3 ↔ 2 process.
To obtain an expression analogous to Eq. (40) under the additional assumption that the correct relic abundance is obtained, i.e. that Eq. (10) is satisfied, we need to understand how the freezeout abundance of DM scales with the model parameters analytically in the KINDER regime. In the next few sections, we will review each thermodynamic phase of the KINDER regime, providing where possible an analytic understanding of the KINDER freezeout process.

Kinetic Decoupling and Cannibalization
As we discussed in Sec. III D, kinetic decoupling occurs at the point when both Eqs. (34) and (35) have just been satisfied. We find that kinetic decoupling is usually controlled by A ↔ f f , i.e. the condition Eq. (34) is fulfilled after Eq. (35). Therefore, for the purpose of analytic estimates, we will assume that this is always true; our numerical results show that elastic scattering can become the process controlling kinetic decoupling at m χ ∼ O(GeV) and large α D .
Let us first obtain an analytic estimate of x d , the dimensionless inverse temperature at kinetic decoupling, to see how it depends on the parameters of our model. Using the expression in Eq. (8) with T = T and µ χ = 0, Eq. (34) reads For our benchmark parameters in this regime, the value of x where this condition is met is shown in the right panel of Fig. 4 at the transition between thermodynamic phases A and B. After kinetic decoupling the dark sector temperature T deviates from the SM temperature T , as indicated in Fig. 3, while the 2 → 2 and 3 → 2 processes continue to proceed at rates larger than the Hubble expansion rate. The dark sector enters thermodynamic phase B, with both the 2 ↔ 2 and 3 ↔ 2 processes maintaining chemical equilibrium in the dark sector and forcing the chemical potentials to zero, as discussed in Sec. III D.
During this phase, the dark sector is cannibalistic, undergoing a net conversion of mass to kinetic energy in the dark sector, which then causes the dark sector to heat up.
In the limit where no energy is transferred to the SM, the dark sector entropy s D a 3 is conserved. The dark sector entropy density can be approximated as where in the second line we can neglect ρ A due to its relatively large Boltzmann suppression compared to ρ χ , and we used the fact that P A P χ = n χ T m χ n χ for x 1. Conservation of entropy enforces d s D a 3 /dt = 0, with no processes active between the dark sector and the SM. In this limit, we have µ χ n χ µ A n A since µ χ and µ A are of the same order, and µ χṅχ + µ A ṅ A 0, since the fast dark sector processes are responsible for both setting the chemical potentials and the number density evolution of the dark sector particles. Making use of Eq. (43) and the expression ofṅ χ in Eq. (8), entropy conservation in the dark sector implies the following relation between T and T : In thermodynamic phase B, we have µ χ ≈ µ A ≈ 0 and m χ T , giving: which we can integrate from T d = T d ≡ m χ /x d up to some dark sector temperature T to get We see that the dark sector temperature T is approximately fixed by the temperature of kinetic decoupling T d , with x evolving slowly (logarithmically) with x thereafter. If entropy were perfectly conserved, then the corresponding evolution in n χ would be approximately which would indicate an approximately constant n χ a 3 and Y χ in phase B. While entropy conservation arguments are sufficient to get a crude approximation of the behavior of the dark sector in this phase, the true picture is significantly more complicated; for example, in Fig. 3, while Y χ stops exponentially decreasing in phase B, it is clearly not constant.
In Fig. 4, we see that the dark sector enters thermodynamic phase B after kinetic decoupling occurs at around x d ∼ 15 for our benchmark parameters. After kinetic decoupling, the A ↔ f f is no longer fast enough to keep the dark sector and SM in thermal equilibrium. As a result, the dark sector begins to heat, as shown in the bottom panel of Fig. 3. With the increase in T , however, comes an increase in n A , which also increases the rate at which energy density is transferred by A → f f to the SM. As a result, the energy density transfer from the dark sector to the SM remains relatively large even after kinetic decoupling; this can be seen in Fig. 4, which shows that this rate stays close to the rate of change of the dark sector energy density per χ particle, given approximately by m χṅχ /n χ . Dark sector entropy is thus not quite conserved.
A better analytic understanding for the dark sector evolution thermodynamic phase B can be obtained from the argument above: since T always evolves in such a way as to keep A → f f relatively efficient at transferring energy from the dark sector to the SM, we find that Thermodynamic phase B is characterized by zero chemical potentials for both species, i.e. n χ ≈ n χ,0 (T ), and likewise for n A . Given the expression forṅ χ in Eq. (8), this approximation gives Taking 3x/2x x and Fig . 5 shows the comparison between this analytic temperature evolution and the numerical evolution computed directly from the Boltzmann equations. We see that the analytic result assuming entropy conservation overestimates the temperature somewhat, since it neglects the transfer of energy to the SM, and our modified analytic estimate is in better agreement with the phase B numerical results. As we indicated earlier, a very similar logarithmic evolution of x in a kinetically decoupled dark sector with zero chemical potential has already been found in other dark sector models [17,18,[29][30][31]. However, as discussed above in Sec. IV B 1, in the dark photon model parameter space we are studying, a second stage of cannibalization begins when the universe expands and cools to the point where the 2 ↔ 2 process freezes out.

2 ↔ 2 Freezeout and Continued Cannibalization
The 2 ↔ 2 process freezes out when the χχ → A A rate falls below the Hubble expansion rate, triggering a nonzero chemical potential in the dark sector; this is indicated on the left panel of Fig. 4 by the transition from phase B to C. We will label the temperatures of the SM and dark sector at which 2 ↔ 2 freezeout occurs as T 2 and T 2 respectively, and correspondingly x 2 and x 2 .
To understand the behavior of the dark sector in this phase analytically, we rely on Eq. (21) and drop the contribution from elastic scattering, which is unimportant by the time the dark sector is in thermodynamic phase C. This giveṡ for the χ number density evolution, anḋ for the A number density.
In general, provided that the dark-sector number densities n i (i = χ, A ) are such that they would be in a steady state in the absence of the cosmic expansion, their time derivatives will be parametrically controlled by H and can be approximated as being of order Hn i (the prefactor, of course, being important to the details of the solution). During the two cannibalization stages, when the comoving number density evolution is slow, we furthermore expect the prefactor to be an O(1) number. Therefore, Eq. (51) shows that: The 3 ↔ 2 term on the right-hand side of Eq. (53) also appears in the Boltzmann equation for A shown in Eq. (15); however, since in general n A n χ , we see that In the parameter space of interest for obtaining the correct relic abundance in thermodynamic phase C, we generally have Γ H by the time T ∼ m χ , as well as As we argued above, we expect the right-hand side of Eq. (52) to be on the order of Hn A ; since both terms on the right-hand side are large compared to Hn A , we expect these terms to be comparable in magnitude. Given that the 3 ↔ 2 rate is on the order of Hn χ as shown in Eq. (53), we therefore arrive at the following important approximate relation that is valid in phase C: How well the last approximation in the equation above is satisfied determines the accuracy of our analytic results: in Fig. 4, we see that this approximation is satisfied up to a factor of 3 throughout phase C.
In thermodynamic phase C with a fast 3 ↔ 2 process, recall from Eq. (32) that the chemical potentials of χ and A are related by µ A ≈ 2µ χ . We can therefore rewrite Eq. (56) as At the point of 2 → 2 freezeout, with the dark sector and SM temperatures being T 2 and T 2 respectively, we still have µ χ (T 2 ) = 0, and so we have from which we finally obtain the following approximate relation for µ χ as a function of T and T : To obtain a full, analytic understanding of the dark sector evolution, we now need to determine T as a function of T after 2 → 2 freezeout. We can once again obtain a rough approximation by taking the dark sector entropy to be conserved, in which case Eq. (44) determines the evolution of T as a function of T . In order to get analytic control of the temperature evolution, we can make the approximations T m χ and µ χ m χ ; the latter condition is true early in phase C since the chemical potential starts at zero. Using the expression for µ χ /T derived in Eq. (59), we find After 2 → 2 freezeout, for values of r that are not too close to 2, we typically have 2(2 − r)x 1, and so we may drop the first term in the equation above to find that We may integrate this approximate expression to obtain . (62) which shows that even during thermodynamic phase C with a nonzero chemical potential in the dark sector, the dark sector temperature T still evolves logarithmically with the SM temperature T . After the freezeout of the 2 ↔ 2 process, the 3 ↔ 2 process alone is sufficient to maintain cannibalization of the dark sector, even though a nonzero dark chemical potential µ has developed. This second stage of cannibalization which occurs in the KINDER scenario is an extension of the conventional cannibalization scenario. It is a critical part of the thermal history of KINDER, because it ensures that after 2 ↔ 2 freezeout and before 3 ↔ 2 freezeout, the dark sector temperature T and comoving number density n χ a 3 continue to evolve slowly, as Fig. 3 shows, remaining mostly fixed by their values at kinetic decoupling. We will explore this slow evolution of n χ in more detail in Sec. IV B 4. As before, entropy conservation is not strictly obeyed due to the fact that A → f f remains quite efficient at transferring energy from the dark sector to the SM; a more sophisticated analytic understanding can once again be attained by examining the Boltzmann equations closely. First, with elastic scattering being unimportant, Eq. (22) shows that there is an approximate relationship between the 2 ↔ 2 and 3 ↔ 2 rates that is applicable even after 2 ↔ 2 freezeout: As we argued in Eq. (53), the 3 ↔ 2 rate is comparable to Hn χ , which leads us to conclude that This expression demonstrates that just after the point of 2 ↔ 2 freezeout, defined in Eq. (31), the χχ → A A and A A → χχ rates remain close to each other, until The fact that these rates are close even after 2 ↔ 2 freezeout can be seen in Fig. 4, immediately after the transition between phases B and C. Before the condition in Eq. (65) is satisfied, we must therefore have µ χ ≈ µ A as well, which together with the fast 3 ↔ 2 requirement that µ A ≈ 2µ χ maintains the chemical potential of the dark sector at approximately zero. Moreover, temperature evolution continues to obey the temperature evolution derived in phase B, shown in Eq. (50). Eventually, n χ decreases to a point where Eq. (65) becomes satisfied at some SM temperature T µ and corresponding x µ ≡ m χ /T µ . Above x µ , the previous argument used to obtain Eq. (59) can be used to obtain a similar expression: and the condition shown in Eq. (56) reduces the χ number density evolution to the following compact form: Using the expression forṅ χ found in Eq. (8) as well as the expression for the chemical potential derived in Eq. (59), we obtain If we make the approximation that 3/2 (2 − r)m χ /T , we can integrate this expression to obtain where Compared to the estimate for x in phase C obtained using entropy conservation in Eq. (62), we see that this more sophisticated analytic treatment (i) correctly identifies the delay in the onset of a nonzero chemical potential, and (ii) introduces a correction to the temperature evolution encapsulated by the factor C (C 3.4 for our benchmark parameters). The result of our analytic estimate for the temperature is shown in Fig. 5, and shows reasonable agreement with the fully numerical solution, up till the complete freezeout of the dark sector at x ∼ 200. The agreement between the analytic estimate and the numerical result deteriorates at larger x as the approximation Hn χ ≈ rΓn A becomes poor (we should only expect them to be equal up to an O(1) factor). The result for our improved analytic estimate for the chemical potentials using Eq. (66) is shown in Fig. 6, and shows good agreement with the numerical results.

3 → 2 Freezeout and Relic Abundance
Cannibalization of the dark sector continues until the universe expands and cools to the point at which 3 ↔ 2 annihilations freeze out at temperature T 3 (and corresponding x 3 ). This marks the freezeout of DM χ, at x 3 ∼ 200 for our benchmark KINDER parameter point, as demonstrated in Figs. 3 and 4. After freezeout, the comoving DM abundance Y χ settles to its constant relic value, and the dark sector temperature begins to evolve as T ∝ T 2 , as expected for a completely decoupled nonrelativistic fluid.
Given the analytic estimates derived in the previous sections, we can now obtain an analytic estimate for the number density of DM particles at 3 ↔ 2 freezeout, given by the condition shown in Eq. (30). We use the assumption of dark sector entropy conservation for simplicity, although a similar conclusion can be reached by using the more accurate analytic results described previously.
The number density of DM at freezeout can be written given the chemical potential in Eq. (59), giving However, the approximate expression for the temperature evolution in thermodynamic phase C found in Eq. (62) allows us to rewrite this as Finally, using the expression for the temperature evolution during thermodynamic phase B in Eq. (50), we can rewrite x 2 in terms of x d , the dimensionless inverse temperature at which kinetic decoupling occurs, giving This remarkable expression shows explicitly that the freezeout abundance is mostly controlled by kinetic decoupling, being exponentially sensitive to x d , up to small power law corrections.
Since the temperature of the dark sector evolves logarithmically after kinetic decoupling, we can make the approximation x 3 ≈ x d in Eq. (73). Substituting the resulting expression into Eq. (30), we obtain the following estimate for n χ at freezeout: We are now ready to obtain an analytic estimate for K/N as shown in Eq. (41), when the "classic NFDM" regime transitions into the KINDER regime in the α Dplane, but now with the requirement that the correct relic abundance is achieved by choosing m χ appropriately at each point in this parameter space. At the regime boundary, kinetic decoupling and 3 → 2 freezeout occur at roughly the same time, i.e. For large values above N/W , the contours follow a constant value of 2 α D , the parameter combination that appears in the expression for σv χχ→f f ; this corresponds to the WIMP regime.
Below N/W , the freezeout of the dark sector transitions into the 'classic NFDM' regime, with the dark sector remaining in thermal contact up till the point of freezeout, and with the abundance controlled solely by when the 3 ↔ 2 process freezes out. Consequentlyas previously discussed in Sec. IV A -the correct relic abundance does not depend on and is only determined by the value of α D , leading to vertical contours.
For yet smaller values of , we eventually encounter the NFDM-KINDER boundary K/N . Within the KINDER regime, the dark matter abundance is determined by the kinetic decoupling process; over much of the parameter space this process is controlled by A ↔ f f , which only depends on , leading to roughly horizontal contours of approximately constant . At larger values of m χ , the elastic scattering process (which depends on α D ) becomes more important, and starts to play a bigger role in determining when kinetic decoupling occurs.
V. 1 r 1.5 We will now focus on the behavior of the dark sector when 1 r 1.5. For these values of r, the 2 ↔ 2 process freezes out after the 3 ↔ 2 process, leading to qualitatively different behavior in the dark sector. Solving the full Boltzmann equations given in Eqs. (14)-(16) reveals a rich and complicated picture, with both the 1. Regime I: the "classic forbidden" scenario. is large enough that A ↔ f f is fast, so that n A n A ,0 (T ); furthermore, χf → χf elastic scattering is sufficiently fast to ensure that the dark sector temperature is nearly equal to the SM temperature throughout the freezeout. The dark sector stays in thermodynamic phase A until the 2 ↔ 2 process freezes out, and no dark sector number-changing processes remain. This regime is precisely the limit studied in Ref. [27].
2. Regime II: n A = n A ,0 (T ), slight cooling. At slightly smaller values of , the process A ↔ f f is still fast enough to maintain n A n A ,0 (T ). However, this condition is insufficient to keep the dark sector in thermal contact with the SM, which cools due to the net conversion of kinetic energy in χ particles into rest mass of the heavier A particles through χχ → A A . In regime II, is large enough for the elastic scattering process, χf ↔ χf , to transfer some heat from the SM to the dark sector, slowing the cooling.
3. Regime III: n A = n A ,0 (T ), rapid cooling. Going to still smaller values of , the A ↔ f f process is still fast enough to lock the number density of A to n A ,0 (T ), but χf → χf is too inefficient to transfer any heat from the SM to the dark sector at any point after 3 ↔ 2 freezeout. In this limit, the rate of cooling is independent of , and the dark sector cools in a manner that only depends on α D .

4.
Regime IV: KINDER. For the smallest values of that we consider, kinetic decoupling of the dark sector from the SM occurs while both the 3 → 2 and 2 → 2 processes have rates that are much faster than Hubble. This shares many of the features of the KINDER regime discussed for 1.5 r 2: the dark sector first enters thermodynamic phase B with zero chemical potential and logarithmic evolution of T with respect to T , and then transitions to thermodynamic phase C after 3 ↔ 2 freezeout.
We will first discuss the broad features of how the dark sector temperature evolves for 1 r 1.5, before examining each of these regimes in turn, focusing on getting some analytic intuition for them. All of our results are once again obtained by solving the Boltzmann equations, Eqs. (14)-(16), numerically.

A. Dark Sector Temperature Evolution
In Regime I, the "classic forbidden" DM regime, the temperature evolution of the dark sector is trivially given by T = T . For the other regimes, the 3 ↔ 2 freezeout divides the dark sector temperature evolution into two important phases.

Temperature Evolution Before 3 ↔ 2 Freezeout
In Regimes II and III, while both the 2 ↔ 2 and 3 ↔ 2 processes are fast, the simultaneous conditions imposed on the chemical potentials shown in Eqs. (32) and (33) are satisfied only if µ χ ≈ µ A ≈ 0. At the same time, is large enough such that n A = n A ,0 (T ); therefore, we must have n χ = n χ,0 (T ) as well, i.e. T = T . Prior to 3 ↔ 2 freezeout, Regimes II and III thus stay in thermodynamic phase A.
For the KINDER-like Regime IV during this phase, the temperature evolution is identical to the KINDER regime with 1.5 r 2, with T = T prior to kinetic decoupling, and the dark sector entering thermodynamic phase B once decoupling occurs. While in thermodynamic phase B, the dark sector particles have zero chemical potential, and the temperature evolves as in Eq. (50).

Temperature Evolution After 3 ↔ 2 Freezeout
Once the 3 ↔ 2 process freezes out, the only process which depletes χ particles is χχ → A A . This process converts lighter χ particles into heavier A particles, removing kinetic energy from the dark sector, resulting in a cooling of the dark sector. The 2 ↔ 2 process enforces µ χ ≈ µ A , which start to take on nonzero values.
As we derived in Sec. III C, the Boltzmann equations enforce certain relations between the rates of the 3 ↔ 2 process, the 2 ↔ 2 process, A ↔ f f and elastic scattering in the nonrelativistic limit. As shown in Eq. (21), we can approximately express the number density evolution of χ particles purely in terms of the elastic scattering rate and the 3 ↔ 2 rate. In Regime III, the number density evolution between 3 ↔ 2 and 2 ↔ 2 freezeout is dominated solely by the 3 ↔ 2 rate, with the elastic scattering term being negligible. Since the 3 → 2 rate has dropped below the Hubble rate in this phase, Regime III is characterized by n χ a 3 being approximately constant, with the dark sector temperature being dependent only on the 3 → 2 rate. In Regime II, the number density evolution is instead dominated by the elastic scattering rate before 2 ↔ 2 freezeout, leading to more rapid evolution of n χ , and less deviation of T from the SM temperature. In the limit of large elastic scattering, n χ → n χ,0 (T ) with T → T , which is the condition found in Regime I.
To understand the behavior of Regimes II and III more quantitatively, we can expandṅ χ in Eq. (21) using Eq. (8) to obtain In Regimes II and III, approximations for µ χ /T after 3 → 2 freezeout can be found. In these regimes, the value of is large enough such that We emphasize, however, that the dark sector temperature T is not equal to T ; rather, the chemical potential µ A evolves in such a way as to maintain the relation above. The χχ → A A process removes kinetic energy from the dark sector, and the exact evolution of T depends on the efficiency of the heat exchange processes between the dark sector and the SM. Writing out the full expression for n A in Eq. (3) and making use of the fact that while the 2 → 2 process is the only process that is fast, Eq. (33) must hold i.e. µ χ ≈ µ A , we find that the chemical potential must satisfy the following relation: Furthermore, the ratio of n χ and n A is completely specified by x since the chemical potentials cancel out: Eqs. (79) and (80) show that given T as a function of T , we will be able to obtain n χ and n A as a function of the SM temperature in Regimes II and III. Eq. (79) provides an expression for µ χ /T , which combined with Eq. (77) gives an expression for T as a function of T after the freezeout of the 3 ↔ 2 process, with n A (T ) ≈ n A ,0 (T ): If we make the further approximation that m χ T, T , this equation takes a particularly simple form, We note that the second term on the right-hand side is typically smaller than the term before it since T m χ , but has been included to improve the accuracy of this analytic result. In terms of x and x, we have The relative importance of each term on the righthand side of Eq. (82), which governs the temperature evolution after 3 ↔ 2 freezeout, separates Regimes I-III. Since the 3 ↔ 2 term is typically less than O(1), the different regimes are distinguished by how large the elastic scattering term is compared to r/(r − 1). In Regime I, throughout the period between 3 ↔ 2 freezeout and 2 ↔ 2 freezeout, we have keeping in mind that n f σvδE χf →χf ∝ (T − T ) (see Eq. (B16) for an expression for σvδE χf →χf ). The fast elastic scattering enforces T T , the assumption of the "classic forbidden" regime. In Regime II, we have instead at some point between the two dark sector freezeout events. In this regime, since n f σvδE χf →χf ∝ (T − T ), the dark sector begins to cool immediately after 3 → 2 freezeout, but once T starts differing significantly from T , the elastic scattering term becomes large enough to slow the cooling process. Finally, in Regime III, between the 3 ↔ 2 and 2 ↔ 2 freezeout events, we always have This is the limit where the elastic scattering process is too inefficient to transfer heat between the two sectors, and therefore the dark sector cooling is rapid and becomes independent of .

B. Regime Boundaries and Characteristics
We will now describe some general characteristics of each regime, providing where we can an analytic description of the dark sector freezeout process. We also explain how to numerically estimate the value of on the α Dplane at which the boundary between the regimes is located.

For
10 −4 , freezeout of the dark sector is controlled by χχ ↔ f f , corresponding to the conventional WIMP regime. For values of smaller than this, we enter regime I, the "classic forbidden" regime, with T ≈ T until the final freezeout of the dark sector. This regime was studied in Ref. [27], where they showed that the dark sector freezeout is determined entirely by when the 2 ↔ 2 freezeout occurs, a purely dark sector process which is independent of .
The "classic forbidden"-WIMP boundary occurs when the 2 ↔ 2 dark sector process freezes out and approximately the same time as χχ → f f , i.e. n χ σv χχ→A A ≈ H ≈ n χ σv χχ→f f (WIMP/I) . (87) If we further require the freezeout to produce the observed relic abundance and fulfil Eq. (10), we obtain the following analytic estimate for WIMP/I , the value of at the WIMP/Regime I boundary, and specializing to r = 1.4 for illustration: where x f ∼ 20 gives the temperature of freezeout of both the 2 ↔ 2 and the χχ ↔ f f processes. At the low-end of Regime I, the elastic scattering energy transfer rate becomes gradually small enough such that Eq. (84) is no longer satisfied at all points between 3 ↔ 2 freezeout and 2 ↔ 2 freezeout, and the dark sector transitions into Regime II. The boundary between Regimes I and II is therefore marked by when the elastic scattering condition for Regime II, Eq. (85), becomes fulfilled just as 2 → 2 freezeout occurs, i.e.
where T 2 and T 2 are the SM and dark sector temperatures at 2 → 2 freezeout. Note that both σvδE χf →χf and σv χχ→A A depend on T and T . Together with Eq. (10) for the relic abundance, we can obtain a numerical estimate for I/II , the value of as a function of α D at the boundary between Regimes I and II.

Regime II
Regime II is characterized by Eq. (85) between 3 ↔ 2 and 2 ↔ 2 freezeout, which ensures that T < T due to the χχ → A A process, but with some heat being transferred from the SM to the dark sector to impede the cooling of the dark sector due to χf → χf . At the same time, the decay rate Γ is large enough such that n A ≈ n A ,0 (T ) throughout the freezeout of the dark sector. This condition immediately determines the chemical potentials µ χ , given analytically by the expression Eq. (79), as well as µ A ≈ µ χ . In Fig. 8, we show this analytic result in comparison with the numeric calculation of the chemical potential, for our Regime II benchmark point of m χ = 10 MeV, r = 1.4, = 10 −6 and α D = 0.03. Note that the agreement deteriorates rapidly once 2 ↔ 2 freezeout occurs at x ∼ 21, after which the DM particle has completely frozen out, and the assumption that µ χ ≈ µ A breaks. A similar result is obtained in Regime III as well, where n A ≈ n A ,0 (T ) also holds. Fig. 9 shows the evolution of the χ number density and T at the same benchmark parameters. n A evolves trivially as n A ,0 (T ), and therefore need not be separately plotted. Since a chemical potential develops immediately after the dark sector kinetically decouples from the SM at the point of 3 ↔ 2 freezeout at x ∼ 15, the dark sector passes from thermodynamic phase A to C directly. The characteristic cooling of the dark sector is apparent in the right panel of Fig. 9, and is governed by Eq. (83). In this regime, this differential equation does not appear to be analytically integrable; we show only the numerical result, obtained directly from the full Boltzmann equations.
In Fig. 10, we show the number density and energy density rates of all relevant dark sector processes. The transition between phases A and C occurs at roughly x ∼ 15, when the backward and forward 3 ↔ 2 rates cease to be approximately equal. This occurs when the 3 → 2 rate is still much larger than the Hubble rate, due to the relation between the rates of the 3 ↔ 2 and 2 ↔ 2 processes enforced by Eq. (22), where the 2 ↔ 2 rate being of order Hn χ allows the 3 ↔ 2 total rate to be much larger than the Hubble rate. Once the dark sector transitions into phase C, we see that the elastic scattering energy density rate per χ particle becomes just a factor of a few smaller than −m χṅχ /n χ , meeting the Regime II criterion laid out in Eq. (85). This shows that a significant amount of heat is transferred from the SM to the dark sector, slowing the cooling rate compared to what happens in Regime III, which we will discuss next.
Within Regime II, as decreases still further, Eq. (85) is met increasingly earlier, leading to a colder dark sector due to the diminishing ability of χf → χf to heat the dark sector. Eventually, the condition Eq. (85) is only met at the point of 3 → 2 freezeout, and no significant amount of heat is transferred to the dark sector after that. This marks the boundary between Regime II and III, i.e.
where T 3 and T 3 are the SM and dark sector temperatures at 3 ↔ 2 freezeout respectively. An analytic estimate for II/III , the value of when these two conditions are satisfied, is Once again, combining the boundary conditions shown above with the observed relic abundance in Eq. (10) allows us to eliminate m χ and x f from the expression above numerically. This numerical expression for II/III forms the boundary between Regimes II and III. Fig. 11 shows the evolution of the χ-abundance and the dark sector temperature in Regime III, for our bench-   Fig. 12, we show the number density and energy density rates per χ particle through the dark sector freezeout. In this regime, the dark sector temperature once again cools rapidly after 3 → 2 freezeout and enters thermodynamic phase C; un-like Regime II, however, elastic scattering plays no significant role in influencing this evolution between 3 ↔ 2 freezeout and 2 ↔ 2 freezeout, as can be seen in the right panel of Fig. 12. The dark sector temperature evolution after 3 ↔ 2 freezeout can be obtained by setting σvδE → 0 in Eq. (83) and neglecting the 2 → 3 rate  (which is much smaller than the forward rate after 3 → 2 freezeout), i.e.

Regime III
Given the approximation for the chemical potential µ χ in Eq. (79), this differential equation can be integrated exactly, starting from x 3 = x 3 , to give In the right panel of Fig. 11, we show this analytic result in comparison with the numeric result obtained from the full Boltzmann equation, and find excellent agreement between them, up to 2 ↔ 2 freezeout at x ∼ 23. Throughout Regime III, n A ≈ n A ,0 (T ) due to the highly efficient A ↔ f f process; as decreases, however, A ↔ f f becomes less and less rapid, and eventually this process becomes too inefficient to keep the dark sector in thermal equilibrium at the point of 3 ↔ 2 freezeout. Below this point, kinetic decoupling between the two sectors occurs before either of the dark sector processes freezes out, leading to the KINDER-like Regime IV. We can estimate the boundary between Regimes III and IV by requiring the 3 → 2 freezeout and kinetic decoupling to occur at the same time, i.e.
These conditions are however identical to the conditions used for estimating the boundary between the KINDER and the NFDM regime for 1.5 r 2 in Eq. (41). This equation can be restated as III/IV ∼ 10 −8 e 9.6(r−1.4) α D 1.0 where III/IV is the value of between Regimes III and IV as a function of various model parameters. Finally, we may once again combine Eq. (95) with the condition for the observed relic abundance in Eq. (10) to numerically derive the boundary between these regimes. Fig. 13 shows the evolution of the χ-abundance and the dark sector temperature in Regime IV, for our benchmark parameters in this regime, m χ = 10 MeV, r = 1.4, = 2 × 10 −9 and α D = 0.6. In Fig. 14, we show the number density and energy density rates per χ particle throughout dark sector freezeout. In Regime IV, kinetic decoupling occurs before either of the 2 ↔ 2 or 3 ↔ 2 processes become slow. This regime is similar to the KINDER regime with 1.5 r 2, exhibiting heating in the dark sector, with the key difference being that the 3 → 2 process is now slower than the 2 → 2 process. In thermodynamic phase A and B, the physics in this regime is identical to that of the KINDER regime with 1.5 r 2, as discussed in Sec. IV B 2. Kinetic decoupling occurs first at a temperature given approximately by Eq. (42), after which the dark sector enters phase B. An approximation for the evolution of T can be obtained by assuming dark sector entropy conservation, leading to

Regime IV
while a more detailed examination of the Boltzmann equations leads to the improved approximation in Eq. (50), i.e.
x ≈ x d Once the 3 ↔ 2 process freezes out, the dark sector enters thermodynamic phase C. As before, the χ number density evolution is given by Eq. For 1 r 1.5, however, the 3 → 2 process is slow in phase C, meaning thaṫ i.e. n χ ∝ a −3 in phase C, with χ frozen out. More accurately, Eqs. (54) and (55) are still true in this regime, since n A n χ ; we therefore still have the following approximate relation after 3 ↔ 2 freezeout occurs: where we have neglected the 2 → 3 rate since the 3 ↔ 2 freezeout has occurred. This approximate relation gives us an expression for µ χ µ A : A comparison between this analytic approximation and the numerical result in phase C is shown in Fig. 15, We can substitute our analytic expression for µ χ /T into Eq. (98) using the expression forṅ χ in Eq. (8), giving where T 3 is the temperature at 3 ↔ 2 freezeout. This expression can be integrated exactly to give This analytic prediction in comparison with the numerical temperature evolution is shown in Fig. 16, showing good agreement until near the 2 ↔ 2 freezeout, when µ χ and µ A begin to diverge.    can be made out by changes in behavior of the contour lines. Note that the boundary between the WIMP regime and Regime I occurs at values above the maximum shown in Fig. 17.
In Regime I, the relic abundance is controlled entirely by the 2 ↔ 2 freezeout, which only depends on α D , leading to vertical contours in the α D -plane. Decreasing into Regime II, the relic abundance is controlled by when the freezeout of 3 ↔ 2 and of 2 ↔ 2 occur, as well as how efficiently χf ↔ χf heats the dark sector and impedes the cooling due to χχ → A A , leading to some nontrivial dependence on and α D . Once we arrive at Regime III however, elastic scattering becomes extremely inefficient, and the rate of dark sector cooling after 3 ↔ 2 freezeout depends only on the 3 → 2 rate itself. Since all of the physically important processes are purely dark sector processes, the contours are once again independent of . Finally, in Regime IV, the relic abundance is determined by when kinetic decoupling occurs, but also by the long power-law decrease in n χ in phase B, which is dictated by dark-sector-only processes. This once again leads to contours that depend on both α D and .
We note that the contour of Ω χ h 2 = 0.12 for m χ 5 MeV shows an abrupt change in behavior in Regime II compared to higher DM masses. This occurs due to the fact that in Regime II thermodynamic phase C, DM particles with masses below ∼ 5 MeV undergo elastic scattering with nonrelativistic, rather than relativistic, electrons throughout most of the freezeout process. The Boltzmann suppression of nonrelativistic electrons leads to a sharp decrease in σvδE , which controls the cooling rate of the dark sector in this phase, and hence the relic abundance of DM. The correct relic abundance is thus achieved at a higher value of than expected, in order for the stronger mixing to compensate for the decrease in electron number density. We refer the reader to App. B for more details on how σvδE χf →χf is computed.

VI. EXPERIMENTAL PROBES AND CONSTRAINTS
There are significant constraints on dark photons from both terrestrial experiments and supernova observations. There are also cosmological constraints on the DM itself, from DM annihilation to electrons and positrons affecting the anisotropies of the CMB, and from modifications to the number of effective degrees of freedom during Big Bang nucleosynthesis (BBN) and the CMB epoch. DM self-interactions mediated by the dark photon exchange can be large, and can be probed by observations of galactic structure. Finally, a sufficiently warm dark sector can be constrained by measurements of the matter power spectrum. We will discuss these constraints in this section, and plot the results in Fig. 18.

A. Accelerator and Direct-Detection Experiments
For 1 r 2 with a dark photon mass above 1 MeV, dark photons produced at beam experiments decay visibly into SM particles. The observational signatures of visibly decaying dark photons have been studied extensively in the literature [6][7][8][9][10][11][12][13][14][15]. In Fig. 18, we plot the region of parameter space excluded by these experiments. This excluded region covers considerable parameter space, extending down to ∼ 10 −7 − 10 −8 for Direct-detection experiments can probe the scattering of the DM on both electrons and nucleons (including the Migdal effect [33,34]) via dark photon exchange. In Fig. 18, we consider the constraints from DarkSide, Xenon 1T, SuperCDMS, and SENSEI [33][34][35][36][37][38]. In the parameter space we consider, nuclear scattering limits derived by exploiting the Migdal effect set the strongest bound. These limits are primarily sensitive to the highmass, high-corner of our parameter space.

B. Supernova Constraints
The production and escape of dark sector particles during a core-collapse supernova can lead to cooling of the proto-neutron star that differs from the SM prediction [39,40]. Such anomalous cooling is constrained by our observation of SN1987A [41,42]. 1 Ref. [44] carefully derived constraints on the m χplane in the vector-portal DM model using the SN1987A result for m χ = 3m A , and for two discrete α D values, together with constraints for models with only A and no DM. For fixed α D , m χ and m A , the excluded region is generally enclosed by two boundary values of . The lower boundary in is determined by the rate of production of the dark-sector particles from the SN core: models with smaller values of are allowed because they do not lead to enough production of dark-sector particles to modify the supernova evolution significantly. The upper boundary on is determined by whether the darksector particles will thermalize with the SM material in the proto-neutron star before escaping the SN, leaving these particles trapped; in models with larger values of , the dark-sector particles are thermalized efficiently and do not escape and cool the proto-neutron star, and hence these scenarios are unconstrained.
We now discuss how to recast the bounds in Ref. [44] for different values of α D . The maximum value of m χ is independent of α D , being set by the kinematics of the supernova. The behavior of the lower bound is determined by the DM mass with respect to the plasma frequency of the interior, ω p ∼ 15 MeV. For 2m χ > ω p , the off-shell DM production via bremsstrahlung through virtual dark photons during neutron-proton collisions is suppressed, and the direct production of A is more important. Consequently, the lower bound in is very similar to that in the dark-photon-only case, and is roughly independent of α D . For 2m χ < ω p , however, χχ-pairs can be produced through an on-shell A , and the production rate is fixed by the value of α D 2 . For a lower bound given at a reference value α D,ref , we can therefore rescale to a new value of α D by leaving the part of the bound where 2m χ > ω p constant, and rescaling the limit where 2m χ < ω p by The upper boundary of the limit on is determined by the dark-matter-proton scattering cross-section, and consequently varying α D changes the asymptotically flat part of the upper boundary in such that α D 2 is kept fixed, i.e. from a reference upper limit given for α D,ref , we rescale by α D,ref /α D .
We find that for 10 −9 , the DM rate of production in the supernova in our model is always large enough for a significant amount to be produced; our limits are therefore set by the upper limit on , as determined by the 1 Alternative cooling models have also been proposed that cast doubt on the SN1987A bounds (see, e.g., Ref. [43]). thermalization condition. Note that this also happens for the lower boundary of our curves since there α D is very large. The SN1987A constraints cover the low-and low-m χ part of the parameter space, and generally lie entirely within the self-interaction constraints that we will describe next (albeit with different model-dependence).

C. DM Self-Interactions
The cross section for elastic DM-DM scattering is constrained by cluster mergers and halo shapes to satisfy σ SI /m χ ≤ 1 cm 2 g −1 ∼ 5 × 10 3 GeV −3 [45]. The DM selfinteraction rates for χχ → χχ and χχ → χχ (and their conjugate processes) are determined in Refs. [27,28]. Including both s and t-channel tree level diagrams, the averaged cross section σ SI is given by: where As shown in Fig. 18, this constraint rules out a large fraction of the parameter space especially at low , generically excluding as high as 10 −6 -10 −5 depending on r; this behavior occurs because the values of α D required to obtain the correct relic density are higher at small . In this sense the self-interaction bound is complementary to limits on the interactions with the SM, which are suppressed by small .
One possible way to evade this constraint is to consider a scenario where only some subdominant fraction of the DM is produced by the mechanisms we have considered in this work, as this limit is rather sensitive to the fraction of DM that is self-interacting. For example, Ref. [46] shows that if the self-interacting component is less than 1% of the DM, these constraints become inapplicable. However, a full self-consistent treatment of fractionally abundance self-interacting dark matter constraints would require recalculation of the cosmological evolution in order to obtain a lower relic density, and is beyond the scope of this work.

D. CMB Constraints on DM Annihilation
During the post-recombination epoch, DM annihilation to e + e − leads to energy deposition into the baryonic gas; the resulting extra ionization can be constrained based on observations of the CMB anisotropy. We compare the annihilation cross section for χχ → ff (see App. B) to the limits derived in Ref. [47] and updated in Ref. [23]. We plot the region excluded by this constraint in Fig. 18.
We observe that these CMB constraints provide some of the strongest bounds on models of this type for r close to 2, excluding most of the available parameter space. Even for smaller values of r, the CMB constraints provide stringent limits for models with low m χ and high .
These limits could be lifted or relaxed if the dark-sector model were adjusted in order to suppress the DM annihilation to SM particles at low velocities. For example, this could be achieved if the DM was a scalar rather than a fermion, as then the leading-order annihilation through the dark photon would be p-wave and scale as σv ∝ v 2 .

E. Cosmological Constraints on Light Relics
Electromagnetically coupled DM with a mass of around 1 MeV can significantly affect the process of Big Bang Nucleosynthesis (BBN) by (i) directly increasing the expansion rate as a contribution to the energy density of the universe, and (ii) injecting entropy into the SM sector and changing the relative energy density of the electromagnetic sector as compared to the neutrino sector, altering the temperature evolution of both sectors with respect to standard cosmology. These changes in turn alter the predicted abundance of light nuclei like deuterium and helium-4, which can then be compared with existing measurements of the abundances of these nuclei (see e.g. Refs. [48][49][50][51][52] for deuterium and helium-4). The injection of entropy from electromagnetically coupled DM can also decrease N eff [53], the effective number of degrees of freedom, during the CMB epoch, which can then be constrained by the CMB anisotropy power spectrum [23].
Ref. [54] modelled the predicted primordial elemental abundances in the presence of an electromagnetically coupled dark matter particle; we adopt their results for our BBN constraints. They presented two constraints, depending on whether a prior was imposed on Ω b h 2 in the BBN calculations. When no prior was imposed, the bound is relatively weak, m χ 0.7 MeV for Dirac fermion DM. With a prior based on CMB observations, Ω b h 2 = 0.02225 ± 0.00066 [23], this bound improves to m χ 7 MeV, since the effect of entropy injection into the SM from the DM cannot be compensated for by lowering Ω b h 2 arbitrarily.
We note however that assuming the central value of Ω b h 2 from Planck leads to a standard BBN theoretical prediction of D/H that is roughly 2σ below the central measured value. This discrepancy may indicate an incomplete understanding of the process of BBN even in standard cosmology, which may therefore affect the bound given above.
As mentioned above, one can also consider the impact of electromagnetically coupled DM particles on the CMB anisotropy power spectrum. Electromagnetically coupled DM particles heat the electromagnetic sector as they become nonrelativistic, effectively decreasing the number of relativistic degrees of freedom at late times by increasing the ratio of photon to neutrino temperatures. The Planck 2018 measurement [23] sets a constraint on electrophilic Dirac fermions of m χ 7.4 MeV. A joint constraint using both primordial elemental abundance and CMB data strengthens the constraint on electrophilic Dirac fermion to m χ 10 MeV. However, CMB N eff bounds are less robust than the BBN constraint, and can be overcome by e.g. adding dark, relativistic degrees of freedom to compensate for the effect of the electromagnetically coupled DM [55].
Given the above consideration, we set a tentative constraint of m χ > 7 MeV to indicate the potential constraint from BBN and CMB. Since the region with m χ < 10 MeV is already strongly constrained by the CMB limits on DM s-wave annihilation, beam dump experiments and SN1987A, this constraint is not particularly important to understanding the viability of the model.

F. Warm Dark Matter
In the 1.5 r 2 KINDER regime, the dark sector undergoes an early kinetic decoupling from the SM, after which the dark sector temperature T evolves only logarithmically with respect to the SM temperature T until the 3 → 2 process freezes out. As a result, the dark sector temperature can be much higher than in the standard WIMP paradigm, where T = T until kinetic decoupling, after which T ∝ (1 + z) 2 . Models of warm dark matter (WDM) typically have suppressed structure on small scales [56,57], and can be constrained by measurements of the matter power spectrum from the Lyman-α forest [58,59], which are sensitive to modes with comoving wavenumber as large as k max ∼ 3 h Mpc −1 .
To get an estimate for how important the WDM Lyman-α bounds are to the KINDER regime, we estimate the comoving Jeans length λ J of DM, and compare this with 2π/k max ∼ 2 h −1 Mpc; for model parameters where λ J 2π/k max , the model is unlikely to leave a significant imprint on the matter power spectrum on scales currently probed by experiments. We leave a detailed analysis of such potential WDM constraints for future work.
The comoving Jeans length for the DM is given by [60] λ J (z) = (1 + z) After the dark sector completely freezes out, T ∝ (1 + z) 2 ; in the radiation dominated era, λ J stays roughly constant, while λ J ∝ (1 + z) 1/2 during matter domination, decreasing with time. To make a conservative estimate, we therefore want to compare λ J (z eq ) with 2π/k max at the redshift of matter-radiation equality, z eq . 2 We can estimate the temperature of the dark sector at z eq as where z 3 and T 3 are the redshift and dark sector temperature at 3 → 2 freezeout respectively. With this approximation, we have Taking z eq = 3402 and assuming a ΛCDM cosmology, we can obtain the following estimate for the Jeans length at matter-radiation equality: In the 1.5 r 2 KINDER regime, we know that x 3 ∼ x d , since x evolves logarithmically with respect to x in thermodynamic phases B and C, while x 3 is largest when the 3 ↔ 2 process freezes out at the latest possible time.
We therefore find that λ J (z eq ) is largest at (i) small , so that decoupling occurs early, minimizing x d and thus x 3 , and (ii) large α D with small m χ , so that the 3 ↔ 2 cross section is large, and the process freezes out as late as possible, maximizing x 3 . To maximize the impact on small-scale structure, we therefore take the smallest mass we consider m χ = 1 MeV, choose the largest perturbative value of α D = 4π, giving = 3.5 × 10 −9 to achieve the observed relic abundance for r = 1.8. We find that x 3 = 5500 and x 3 = 45, leading to λ J (z eq ) 0.5 h −1 Mpc, which is still small enough to be consistent with probes of small-scale structure. Other parameter combinations that obtain the observed relic abundance lead to smaller values of λ J (z eq ).
For 1 r 1.5, Regime I has T = T until freezeout of the dark sector, while in Regimes II and III, the dark sector is actually colder than a dark sector that is thermally coupled to the SM until freezeout, easily avoiding these warm DM constraints. In Regime IV, a similar argument as above shows that λ J (z eq ) is given by Eq. (109) with x 3 , x 3 replaced by x 2 , x 2 . Once again, large values of α D , small values of and small m χ would lead to the largest impact on small-scale structure. SN1987A constraints and the requirement of a perturbative value of α D < 4π, however, are enough to constrain 10 −9 . Choosing r = 1.4, = 10 −9 , m χ = 1 MeV and α D = 0.19, we find x 2 = 630, x 2 = 89 and λ J (z eq ) 5 × 10 −3 h −1 Mpc, much smaller than would be observable. Larger values of m χ require larger values of to meet the relic abundance criterion, and lead to even smaller values of λ J (z eq ). Similar results hold for r = 1.3 as well.
We therefore find that λ J (z eq ) 2π/k max is satisfied throughout all relevant parameter space, leaving our model unconstrained by small-scale structure observations. However, parts of the KINDER regime are close to being constrained by existing power spectrum measurements; future improvements in WDM constraints could potentially probe these models. Fig. 18 shows a plot of the constraints on the m χplane with four different values of r, with α D chosen at every point in parameter space such that the observed relic abundance of DM is attained, Ω χ h 2 = 0.12. Regions ruled out by the constraints discussed above are marked in color; parts of the space that require α D > 4π to obtain the correct relic abundance are also shaded gray, since perturbative control of our model breaks down there. The contour of α D = 1 is also shown for reference. For 1.5 r 2, we show the constraints for two representative values, r = 1.6 and r = 1.8. In both cases, a small region of open parameter space exists near ∼ 10 −6 and with DM masses of a few hundred MeV. For these values of r, the vector-portal DM model is bounded from below by the nonperturbative region, and is strongly constrained by the CMB s-wave annihilation bound and selfinteraction limits. The available parameter space sits in the NFDM regime for r = 1.6, and in the KINDER regime for r = 1.8. The unconstrained regions are similar to those obtained in Ref. [28] at the high-end, but differ at the low-end due to the KINDER regime that we have found in this paper.

G. Summary of Constraints
For 1 r 1.5, we show the constraints for r = 1.3 and r = 1.4. Here, there are two viable regions of parameter space: both are in the range m χ 100 MeV, and are separated by the beam dump constraints: one region in Regime III is in the range ∼ 10 −8 -10 −7 , while the other is in Regime II and I in the range ∼ 10 −7 -5 × 10 −5 . In this range of r-values, both the selfinteraction and CMB s-wave annihilation limits are less constraining, allowing more open parameter space than for 1.5 r 2. These new limits represent an improved calculation over those found in Ref. [21]. In particular, most of the available parameter space is not in Regime I, as assumed by Ref. [21]. In contrast to that work, we find that there is a lower limit of 10 −8 imposed by perturbativity and self-interaction constraints, since (in Regime IV) α D needs to become very large at such small values of in order to achieve the correct DM relic abundance.
We emphasize that these constraints are derived assuming that the dark sector is in thermal equilibrium with the SM at T ∼ m χ , which may not be a valid assumption for values smaller than eq as defined in Eq. (37). For ∼ 10 −9 and below, other mechanisms such as freeze-in can potentially achieve the correct relic abundance without the dark sector ever being in thermal equilibrium with the SM.
H. Lifting CMB and Self-Interaction Constraints with Pseudo-Dirac DM In the previous subsections, we have demonstrated that the bulk of the parameter space for this class of models with 1.2 < r < 1.8 has been tested by existing observations and experiments, for the baseline scenario where the DM is a Dirac fermion. Narrow regions of parameter space remain open, but for example, Regime IV for 1 r 1.5 appears to be fully excluded.
However, these exclusions rely critically on constraints from the CMB and from self-interactions, both of which probe the behavior of the DM long after freezeout. This exclusion is model-dependent; it is possible to perturb our baseline model in ways that dramatically alleviate these constraints while leaving the cosmology during the freezeout epoch essentially unchanged.
As a specific example, suppose that the DM is a pseudo-Dirac fermion, where at low energies the DM is split into two nearly-degenerate Majorana mass eigenstates χ 1 , χ 2 (see e.g. Refs. [61][62][63] for specific models). The gauge interaction between the DM and the A (χ / A χ) then gives rise to interactions of the formχ i / A χ j , i = j. There is noχ i / A χ i vertex as Majorana fermions cannot carry a conserved dark charge. The heavier mass eigenstate χ 2 can thus decay to the lighter eigenstate χ 1 via emission of an off-shell A .
When the temperature of the dark sector exceeds the mass splitting between the states, the DM will behave as a Dirac fermion, and thus for a mass splitting ∆m χ T throughout freezeout, our previous cosmological results will still hold. However, once T ∆m χ , the DM will convert into the lighter mass eigenstate provided the lifetime of the heavier eigenstate is sufficiently short (even if the lifetime is long, DM-DM scattering can also efficiently deplete the heavier eigenstate). Thus during the recombination epoch and in galaxies at late times, any process requiring the presence of both mass eigenstates will be strongly suppressed.
This suppression applies to both the annihilationχχ → e + e − through an s-channel A , which determines the CMB constraint, 3 and to the contribution to the treelevel self-interaction cross sectionχχ →χχ from an s-channel A . The contribution to the tree-level selfinteraction cross sections from a t-channel A exchange is suppressed for a related reason; if the initial state is χ 1 χ 1 then the final state (at tree level) can only be χ 2 χ 2 , which is kinematically forbidden provided the kinetic energy of DM particles in the halo is much smaller than the mass splitting. There will still be a contribution to the self-interaction cross section at 1-loop order, and a CMB signal via t-channel annihilation of χ 1 's to the 3-body final state A + e + + e − [64] (as well as possible contributions from the residual χ 2 abundance), but these rates are parametrically suppressed compared to those relevant for the Dirac case.
Thus we expect both the CMB and self-interaction limits to be dramatically relaxed in the pseudo-Dirac case without changing the freezeout history, for mass splittings that are small compared to T at freezeout, but large compared to the DM temperature during recombination and the kinetic energy of DM particles in presentday halos. This modification opens up allowed parameter space spanning all the freezeout regimes we have studied; we will present a detailed computation of the modified constraints in future work [65].

VII. CONCLUSION
We have fully characterized the possible freezeout histories of the vector-portal DM model in Eq. (1), in the region of parameter space in which the DM is a thermal relic, and 1 r 2. In this region, the χχχ ↔ χA (3 ↔ 2) and kinematically suppressed χχ ↔ A A (2 ↔ 2) processes play important roles in the thermal freezeout of the DM. Extending beyond the scope of previous studies [21,22], we explored this model for values of the kinetic mixing parameter where the dark and SM sectors do not remain in kinetic equilibrium throughout the process of DM thermal freezeout. Doing so reveals a rich set of novel thermal histories, leading to very different dependences of the DM relic abundance on the model parameters.
We have identified four novel pathways by which thermal freezeout of the dark sector can proceed, in addition to those identified in previous studies. Two of these pathways share key features, and represent a general class of freezeout histories that we dub the "KINetically DEcoupling Relic" (KINDER). In the KINDER scenario, the DM relic abundance is determined primarily by the kinetic decoupling of the dark and SM sectors. KINDER is realized through a process of dark sector cannibalization, which was previously invoked in the ELDER scenario [17,18]. In this work, we have demonstrated that cannibalization can be supported by a 3 → 2 annihilation process involving multiple dark sector species, and can proceed even in the presence of nonzero dark sector chemical potentials. ELDER DM can be regarded as an example of a KINDER scenario where the kinetic decoupling is controlled by elastic scattering between the DM and SM.
We have presented detailed numerical results for the thermal history of the dark sector in each of these new regimes. Additionally, in a number of cases we were able to analytically derive the evolution of the dark sector temperature T and dark matter abundance Y χ , throughout the freezeout of the DM; this allows us to analytically demonstrate the dependence of the DM relic abundance on the model parameters in much of parameter space.
The novel freezeout mechanisms we have characterized, and their corresponding distinct regimes of parameter space, can be separated into two main parameter regions in r. In the region 1.5 r 2, in addition to the "classic not-forbidden" regime studied in Ref. [22], we have identified a realization of KINDER at low values of .
In the region 1 r 1.5, in addition to the "classic forbidden" regime studied in Ref. [21] (Regime I), which is valid at high , we identify a second variation of KINDER at very low (Regime IV). At intermediate values of , we find two previously unrecognized parameter regimes with distinct freezeout histories (Regimes II and III). In Regimes II and III the A → ff process is fast enough to maintain n A ≈ n A ,0 (T ) until all number-changing processes have frozen out. However, during the period after 3 → 2 freezeout and before 2 → 2 freezeout, this process cannot maintain thermal equilibrium between the DM and SM sectors due to number and energy conservation requirements enforced by the Boltzmann equations. In these regimes the elastic scattering χf → χf process controls the heat exchange between the DM and SM sectors after the freezeout of the 3 ↔ 2 process and before the freezeout of the 2 ↔ 2 process, while the 2 ↔ 2 process cools the dark sector.
The distinguishing feature between Regimes II and III is the efficiency with which the elastic scattering process heats the dark sector. In Regime III, elastic scattering is inefficient, the dark sector is cooled by the kinematically forbidden χχ → A A (2 → 2) process, and the chemical potential of the dark sector is such that the χ abundance no longer evolves appreciably after the 3 ↔ 2 process freezes out. This leads to a DM relic abundance determined only by the freezeout of the 3 ↔ 2 process, even though the 2 ↔ 2 process is significantly faster. In Regime II, in contrast, elastic scattering remains efficient after the freezeout of the 3 ↔ 2 process, and can counteract the cooling of the dark sector, allowing continued evolution of the DM density. This leads to a DM relic abundance determined by the interplay of elastic scattering and dark sector processes.
The two variations of KINDER we have identified differ in their evolution at late times, after the slower dark sector process freezes out. For 1.5 r 2, cannibalization continues through the 3 → 2 process until all number-changing processes have frozen out, ensuring a slow evolution of the DM number density after kinetic decoupling. In contrast, for 1 r 1.5, the cannibalization is halted once the 3 ↔ 2 process freezes out. The number-changing 2 ↔ 2 process is still active at this point, and cools the dark sector; however, the chemi-cal potential evolves such that the χ abundance remains constant regardless.
We have calculated the relevant experimental constraints on our model. Our results drastically modify those of Ref. [21] for 10 −5 (below Regime I) and those of Ref. [22] for 10 −7 (the NFDM and KINDER Regimes). The KINDER mechanism realized in our model implies large self-interaction rates, and a large s-wave annihilation signal in the CMB, for symmetric Dirac fermion DM; these limits are in tension with the KINDER regime, although a small window of open parameter space remains for r = 1.8. There is also available parameter space in Regimes II and III for DM masses ∼ (0.1 − 1) GeV where experiments have not yet explored. In these allowed regions of parameter space, self-interactions can be in the correct range (0.1 cm 2 /g σ SI /m χ 1 cm 2 /g) to have observable consequences for the small-scale structure of galaxies without being currently excluded. Our new calculations provide target regions that can be tested by future sub-GeV direct detection experiments and dark photon searches.
In this paper we have presented the baseline scenario of this vector-portal model in which the DM χ is a Dirac fermion. In a forthcoming paper [65], we will present an alternative to this baseline scenario in which the DM is a pseudo-Dirac fermion which at low energies splits into two nearly-degenerate Majorana mass eigenstates. For the correct range of values of the mass splitting this scenario shares essentially the same cosmology as the Dirac case, while modifying the late-time cosmology in a way that relaxes both CMB and self-interaction constraints, thus opening windows of parameter space spanning all the novel freezeout regimes we have presented.