GeV-scale dark matter with p -wave Breit-Wigner enhanced annihilation

We consider a light scalar dark matter candidate with mass in the GeV range whose p -wave annihilation is enhanced through a Breit-Wigner resonance. The annihilation actually proceeds in the s -channel via a dark photon mediator whose mass is nearly equal to the masses of the incoming particles. We compute the temperature at which kinetic decoupling between dark matter and the primordial plasma occurs and show that including the effect of kinetic decoupling can reduce the dark matter relic density by orders of magnitude. For typical scalar masses ranging from 200 MeV to 5 GeV, we determine the range of values allowed for the dark photon couplings to the scalar and to the standard model particles after requiring the relic density to be in agreement with the value extracted from cosmological observations. We then show that µ and y -distortions of the CMB spectrum and X-ray data from XMM-Newton strongly constrain the model and rule out the region where the dark matter annihilation cross-section is strongly enhanced at small dispersion velocities. Constraints from direct detection searches and from accelerator searches for dark photons offer complementary probes of the model.


I. INTRODUCTION
Searches for dark matter (DM) have for the last decades concentrated on a new weakly interacting particle at the electroweak scale.This was partly motivated by theoretical preferences for new physics at the TeV scale and by the fact that the freeze-out mechanism for a new weakly interacting massive particle (WIMP) at the electroweak scale naturally leads to the value for the DM relic density extracted from observations [1,2].The lack of signals in (in-)direct detection or at colliders for new particles [3][4][5] aroused interest for a wider range of DM candidates of different scales and/or interaction strengths.In particular, on the theoretical side a plethora of DM models have emerged that feature the possibility of lighter DM, around the GeV scale or below [6][7][8][9][10].These candidates largely escape the strongest constraints from direct detection since the minimum recoil energy they can impart to the nucleus falls below the detector threshold [11][12][13].Efforts to improve the sensitivity for direct detection of light DM are being pursued by taking advantage of scattering on electrons or single phonon excitations in crystals [14].Indirect searches for DM annihilation in the galactic halo or in dwarf spheroidal galaxies (dSPhs) are also mostly sensitive to DM masses above a few GeVs [15][16][17].However, WIMPs in the MeV-GeV scale can be constrained from cosmic-ray electrons and positrons using Vogager 1 and AMS-02 data [18,19].Moreover, light DM annihilating into e + e − pairs will leave a signature in X-ray when the electron-positrons scatter on the low-energy photons in the Galaxy.It was shown recently that X-ray data from XMM-Newton can be used to constrain light DM [20,21].
Strong constraints on light DM come from cosmology.DM annihilation deposits electromagnetic energy in the primordial plasma and impacts the anisotropies of the cosmic microwave background (CMB) [22].The precise measurements of the CMB by the PLANCK satellite thus put robust constraints on DM annihilation cross-sections into photons or charged particles.These constraints typically exclude the cross-section required to achieve the measured relic density when the DM mass lies below approximately 10 GeV under the assumption of s-wave DM annihilation [23].However, such constraints are escaped easily if DM annihilation is p-wave, indeed the typical DM velocity during recombination (v ≈ 10 −8 c) is much smaller than the typical velocity at freeze-out, leading to a strong suppression of the cross-section at recombination.The energy injected in the primordial plasma from DM annihilations can also induce deviations from the pure black-body spectrum of the CMB.Depending on the value of the redshift at which the energy injection occurs, the dominant effect is either µ-distortion or, at lower redshifts, y-distortion [24].FIRAS [25,26] has set limits on both types of distortions resulting in important constraints for p-wave annihilating The Lagrangian which we consider describes a vector boson portal, between DM and the SM sector, in the presence of a Breit-Wigner resonance which may enhance DM annihilations.This Lagrangian contains the couplings of a new photon A ′µ , dubbed dark, to ordinary fermions f and to a new complex scalar field ϕ which plays the role of the DM candidate.The scalars ϕ and their anti-partners φ are charged under the new local gauge group U ′ (1) whose associated gauge boson is the dark photon A ′µ .The model to be explored is described by the interaction Lagrangian The DM species ϕ and φ interact weakly on short distances through the exchange of the massive vector boson A ′µ .They respectively carry the dark charges +g x and −g x .The hidden sector has a broken U ′ (1) symmetry, which facilitates the rotation between true A µ and dark A ′µ photon states.This explains the presence of the mixing angle ϵ with which the dark photon couples to ordinary fermions in the second term of the right-hand side of the previous equation.To stabilise the DM, we further impose a discrete Z 2 symmetry, where we assume that ϕ and φ are odd under Z 2 , whereas A ′ µ and all other SM particles are even under the symmetry. 1e can readily extract from the Lagrangian II.1 the cross-section of the annihilation of DM species into ordinary fermions through the reaction The annihilation is mediated by the dark photon A ′µ exchanged in the s-channel.Let us denote by v rel the difference v ϕ − v φ between the ϕ and φ velocities.A straightforward calculation leads to where v is the norm ||v rel || of the velocity difference v rel .The masses of the DM scalars ϕ and φ, and of the dark photon, are respectively denoted by m ϕ and m x , while Γ x is the decay width of the dark photon.The effective dimensionless charge Q is defined by where the sum runs over the fermions f produced by the annihilation, i.e. the fermions whose mass m f is less than √ s/2.The incoming species ϕ and φ are non-relativistic, i.e. v is very small compared to the speed of light.As a consequence, the square of the center-of-mass energy √ s reduces to and, in the non-relativistic limit, the effective charge Q can be defined as For the relic density computation, we are interested in the product σ ann v averaged over the velocity distributions of ϕ and φ species.There is no asymmetry between these components and the velocity distributions for v ϕ and v φ are taken identical.The thermal average of σ ann v is defined as As discussed in Sec.III, we assume hereafter that velocities are distributed according to a Maxwellian distribution where Σ is the one-dimensional dispersion velocity.The integral II.7 can be carried out analytically by swapping the velocities v ϕ and v φ for the variables (II.13) The dark photon decay width Γ x can be expressed as (II. 14) The terms inside brackets respectively refer to decays into ϕ φ pairs and into fermions.The new effective dimensionless charge Q ′ is defined as (II. 15) In our parameter region of interest, m x is slightly larger than 2m ϕ and Q ′ boils down to the charge Q as defined in relation II.6.
From the definitions of a and b, it is obvious that the sign of a determines whether annihilation takes place above or below the resonance.For both cases, however, we can obtain interesting approximations while computing J(a, b) and eventually, ⟨σ ann v⟩.The details regarding this are given in Appendix A.
In this work, we concentrate on the case where the dark photon mass m x is larger than 2m ϕ .Fig. 1 shows how the annihilation cross-section ⟨σ ann v⟩ depends on the one-dimensional dispersion velocity Σ when a < 0, i.e. when the mass gap ∆ = m x − 2m ϕ is positive.As long as ∆ is small with respect to m ϕ , the mass degeneracy parameter Σ 2 0 is nearly equal to the ratio ∆/m ϕ .Substituting Σ from Eq. II.11 into relation II.10, one obtains ⟨σ ann v⟩ = g 2 x ϵ 2 e 2 Q2 3πm 2 x Σ 2 0 × {|a|J(a, b)} , (II.16) In Fig. 1, we plot the term inside brackets as a function of 1/ |a| ≡ (2m ϕ /m x )(Σ/Σ 0 ).For Λ 0 values that are large with respect to Σ 0 , i.e. when Γ x /2 is larger than the mass gap ∆, the cross-section is p-wave suppressed at low velocities and increases like Σ 2 .It reaches its maximum for a velocity of order Λ 0 above which it decreases like Σ −2 .At high velocities, ⟨σ ann v⟩ behaves actually as if it was Sommerfeld enhanced.The model has this intriguing and fascinating property to predict together, albeit in different velocity regimes, a p-wave suppression as well as a Sommerfeld-like behavior.For values of Λ 0 small with respect to Σ 0 , i.e. when Γ x /2 is smaller than the mass gap ∆, a Breit-Wigner resonance comes also into play, which increases dramatically the annihilation cross-section as featured by the bumps of Fig. 1.The maximum is reached at a velocity of Σ M = 2/3 (m x /2m ϕ ) Σ 0 , above which ⟨σ ann v⟩ behaves as Σ −3 .At low and high velocities, the annihilation cross-section scales respectively as Σ 2 (p-wave) and Σ −2 (high-velocity) like in the large Λ 0 case.

III. SCALAR DARK MATTER RELIC ABUNDANCE
From now on, we will assume that scalar DM is in thermal and chemical equilibrium with the primordial plasma when its temperature is of order the DM mass m ϕ .This is the starting point of our analysis.Going beyond this assumption would require the knowledge of the complete theory in order to determine the thermal behavior of DM at much earlier times than those considered here.Because there is potentially an infinite number of such theories, the scope of this article is to focus on the cosmological consequences of a vector boson portal with Breit-Wigner resonance as encoded by the Lagrangian II.1.The model has only a few parameters, namely the dark charge g x , the mixing angle ϵ and the masses m ϕ and m x , but its phenomenology is already very rich and subtle as we shall see.0 /Σ 2 0 .On the horizontal axis, the rescaled variable (2m ϕ /mx)(Σ/Σ0) has been used.When Λ0 is smaller than Σ0, the cross-section is enhanced by a Breit-Wigner resonance.Above a velocity of order Σ0, where its peak value is reached, ⟨σannv⟩ drops like Σ −3 to reach the asymptotic behavior Σ −2 .Below the peak, the p-wave annihilation regime sets in and ⟨σannv⟩ is proportional to Σ 2 .For large values of Λ0 with respect to Σ0, the two asymptotic regimes only appear.
Assuming that scalar DM is in thermodynamical equilibrium with the SM plasma while becoming non-relativistic is quite conceivable.To commence, in the portion of the parameter space where the couplings g x and ϵ are not too small, scalar DM collides upon, and annihilates into, SM fermions so efficiently that thermodynamical equilibrium ensues at early times.If this is not so, we can assume that the kinetic mixing ϵ between the visible and dark photons is triggered by radiative corrections implying loops of heavy dark species Ψ and Ψ ′ with opposite electric or dark charges.At very early times, these particles are relativistic and interact efficiently with both scalar DM and SM species through collisions and annihilations, allowing for the coupling between these components.Notice also that the bulk of DM annihilation takes place well after freeze-out and is most efficient at late times, when the Breit-Wigner resonance becomes active.The DM density at freeze-out turns out not to be relevant to compute the relic abundance.This alleviates the problem of determining the actual thermodynamical state of DM before freeze-out.
We will also assume that scalar DM reaches at all times inner thermal equilibrium so that its velocity distribution is well described by a Maxwell-Boltzmann law.Such a condition could be established through collisions of ϕ and φ particles with SM fermions, provided that energy is transferred sufficiently rapidly between these components.But this turns out not to be the case in a large portion of the parameter space where DM kinetic decoupling occurs.To simplify an already intricate analysis, we have assumed that the ϕ and φ populations reach thermal equilibrium through mutual collisions.This allows us to define a common temperature T ϕ for scalar DM which may be different from the plasma temperature T after kinetic decoupling has occurred.Notice that if the dark charge g x is not too small, DM scalars ϕ and φ could efficiently collide upon each other through the exchange in the t-channel of a virtual dark photon.DM needs to be dense enough though.Alternatively, a self-coupling à la ϕ 4 could also lead to such a behavior.The annihilation cross-section of scalar DM into SM species can then be averaged over a Maxwell-Boltzmann distribution of DM velocities, yielding the result of Sec.II.
In this section, we derive the relic abundance Ω ϕ h 2 of scalar DM as a function of the parameters of the model.
In Sec.III A, we discuss the master equation that drives DM freeze-out and explain how we solve it.We then present numerical results together with approximate solutions that help understanding how Ω ϕ h 2 depends on the mass degeneracy parameter Σ 2 0 .We find three different regimes.In Sec.III C, we assume that scalar DM stays in thermal contact with the SM plasma throughout the entire freeze-out process.This assumption is abandoned in Sec.III D where we determine when scalar DM decouples from kinetic equilibrium and analyse how this strongly affects the relic abundance Ω ϕ h 2 .
A. Calculation of the relic abundance Ω ϕ h 2 The scalar DM particles ϕ and φ can annihilate into, and be produced from, SM fermions through the process Should direct and reverse reactions be fast enough, a chemical equilibrium is established.When the plasma temperature T drops below the DM mass m ϕ , ϕ φ pairs are less and less easily produced.DM still goes on annihilating until it is so depleted that its codensity remains constant until today.We assume that there is no asymmetry between the ϕ and φ populations and denote by n their densities n ϕ = n φ.The evolution of the density n with time t is described by the freeze-out equation where H is the expansion rate of the Universe.In its right-hand side, the three terms respectively stand for the dilution resulting from the expansion of the Universe, the annihilations of ϕ φ pairs and the back reactions which regenerate scalar DM from lighter species.We have assumed detailed balance for this last process, hence the density n eq as fixed by the thermodynamical equilibrium.It is convenient to deal with the codensity ñ, i.e. the density inside a volume that expands with the expanding Universe.We define it as the ratio where θ is a fictitious temperature which decreases as the inverse a −1 of the scale factor a of the Universe, and which is normalized in such a way that it is equal to the plasma temperature T F at freeze-out.Notice that since the entropy of the plasma is constant in time, the factor θ 3 is proportional to the entropy density s = (2π This form allows for a simple interpretation of freeze-out as a mere relaxation process during which ñ runs after its chemical equilibrium ñeq and eventually fails to reach it.The rate Γ F rel ≡ ⟨σ ann v⟩n at which ñ relaxes toward ñeq is actually the annihilation rate.The right-hand side term of relation III.A straightforward calculation yields where the parameter x denotes the DM mass to plasma temperature ratio m ϕ /T .The Hubble expansion rate H is related to the plasma temperature through where M P is the Planck mass and g eff (T ) is the effective number of energetic degrees of freedom of the plasma.In deriving III.6, we have noticed that d ln θ/dt ≡ −H and used the identity The derivative of the annihilation cross-section with respect to the temperature depends on the one-dimensional DM dispersion velocity Σ.In the high-velocity regime where Σ is well above Σ 0 , the cross-section ⟨σ ann v⟩ scales as 1/Σ 2 , hence a value for d ln ⟨σ ann v⟩/d ln x of 1.In the Breit-Wigner enhancement regime just above Σ 0 , this derivative reaches 3/2 while in the low-velocity p-wave regime, it decreases down to −1.
At high temperature, the relaxation rate Γ F rel is much larger than the evolution rate Γ F eq .The chemical equilibrium III.1 is established.The DM density n relaxes toward its equilibrium value n eq much faster than the latter evolves.As time goes on, the relaxation/annihilation rate Γ F rel decreases very rapidly.It actually drops exponentially with plasma temperature insofar as it is proportional to the DM density The rate Γ F eq decreases approximately like T , at a much slower pace than Γ F rel ∝ exp(−m ϕ /T ).At some critical temperature T F , both rates are equal.Freeze-out occurs.Scalar DM decouples from chemical equilibrium insofar as its density n no longer relaxes toward n eq .The freeze-out point x F ≡ m ϕ /T F satisfies the equality This equation must be solved numerically for each set of model parameters, i.e. once the masses m ϕ and m x , the coupling g x and the mixing angle ϵ are given.Expressions II.10, II.11 and II.13 are used to calculate the average annihilation cross-section ⟨σ ann v⟩.Assuming that scalar DM is in thermodynamical equilibrium with the primordial plasma until it becomes non-relativistic imposes that x F is larger than 1.We checked that the solution of equation III.10 actually fulfills this condition in a large portion of parameter space, and noticeably in the domain that survives our cosmological analysis.After freeze-out, the density n eq vanishes rapidly and so does the right-hand side of relation III. 4. We can neglect it in deriving the final DM abundance, especially as most of the annihilation takes place well after freeze-out.The DM codensity today ñ0 is given by the relation The codensity at freeze-out is ñF ≡ ñe (T F ). Equation III.4 is integrated without its right-hand side from time t F , at which freeze-out occurs, until the present age t 0 of the Universe.The annihilation integral I ann may be conveniently recast in terms of an integral over the parameter y = T /m ϕ , where T is the SM plasma temperature.This yields dy ⟨σ ann v⟩M P m ϕ g eff (T ) (III.12) The parameter y runs from y F = 1/x F down to y 0 = T 0 /m ϕ , where the CMB temperature T 0 = 2.72548 K has been borrowed from [40].Once the DM codensity at the present epoch is known, we readily derive the number density and DM mass density ρ 0 ϕ = 2 m ϕ n 0 , remembering that there are as many ϕ as φ particles.Scalar DM contributes today to the Universe mass budget a fraction Newton's constant of gravity is G while H 0 stands for a fiducial Hubble expansion rate of 100 km/s/Mpc.The actual Hubble constant is equal to h in units of that benchmark value.We would like DM to be made of the scalar species ϕ and φ.That is why we will look in Sec.IV for configurations of parameters for which Ω ϕ h 2 is equal to the cosmological measured value of Ω DM h 2 = 0.1200 [1].A key ingredient in the calculation of I ann is the dependence of ⟨σ ann v⟩ on the one-dimensional dispersion velocity Σ of scalar DM, and eventually on the parameter y.Scalar DM is assumed to have reached inner thermalization with temperature T ϕ , hence the identity Σ 2 ≡ T ϕ /m ϕ in the non-relativistic limit.But the relation between DM temperature T ϕ and plasma temperature T remains to be determined.It will prove to be of paramount importance.
B. Approximate expression for the relic abundance Ω ϕ h 2 Before solving numerically for Ω ϕ h 2 , we derive approximate expressions for I ann .These will be helpful to discuss and understand our numerical results.We will set the DM mass m ϕ , the coupling g x and the mixing angle ϵ constant and will concentrate on how I ann , and eventually the DM relic abundance, vary with the mass degeneracy parameter Σ 2 0 ≡ 1 − 4m 2 ϕ /m 2 x .First, we remark that the annihilation cross-section ⟨σ ann v⟩ is not constant after freeze-out.As discussed in Sec.II, it increases with time like 1/Σ 2 in the high-velocity regime, or 1/Σ 3 at the Breit-Wigner resonance, as plasma and DM temperatures decrease.It reaches a maximum when the DM dispersion velocity Σ is equal to some critical value Σ M .This occurs at time t M when the plasma temperature is T M and the parameter y is equal to y M = T M /m ϕ .Afterward, the annihilation becomes p-wave dominated and ⟨σ ann v⟩ drops rapidly like Σ 2 to vanish.In expression III.12, the main contribution to the integral I ann is provided by values of y between y M and y F .From time t M until today, DM essentially does not annihilate.We also remark that the thermodynamical coefficients g eff , h eff and d ln h eff /d ln T vary slowly in time.They can be evaluated at time t M and taken out of the integral over y.Finally, as most of DM annihilation takes place well after freeze-out, the DM codensity ñF is much larger than the present codensity ñ0 and can be removed from expression III.11.Taking into account these remarks and inserting expression II.10 into integral III.12 yields where the plasma function P is defined as The DM relic abundance may be approximated by where It is proportional to the inverse of J ann .Understanding how this integral varies with Σ 2 0 is paramount to our problem.
The ratio b/|a| is the cornerstone of our analysis.It controls how J(a, b) evolves with Σ 2 and, in fine, with y.If that ratio is smaller than 1, a Breit-Wigner enhancement appears with ⟨σ ann v⟩ scaling like 1/Σ 3 for dispersion velocities just above Σ 0 .On the contrary, if it exceeds 1, this enhancement disappears and the annihilation crosssection ⟨σ ann v⟩ is largest for dispersion velocities of order Λ 0 .The ratio b/|a| determines also the point y M where the integrand J(a, b)/Σ 2 in the integral J ann is maximal.It can be expressed in terms of the parameters of the model as For illustration, the ratio b/|a| as a function of Σ 2 0 is presented in Fig. 2 for the same parameters as in Fig. 3.The ratio is minimal when the mass degeneracy parameter Σ 2 0 is equal to the special value This corresponds to a minimal ratio of The mixing angle ϵ cannot be larger than 1 by construction.The dark charge g x is also restricted to be smaller than 1 since otherwise, the theory becomes non-perturbative.According to relation II.15, the effective charge Q ′ 2 reaches a maximal value of 8 in the unrealistic situation where the dark photon mass m x becomes infinite.Even in this extreme case, the minimal value of the ratio b/|a| is less than 1.8 × 10 −2 .As will be discussed in Sec.IV D 2, collider experiments yield furthermore the constraint ϵ ≤ 10 −3 on the mixing angle for m x ≤ 10 GeV [36].This translates into an upper limit on the minimal value of b/|a| of order 1.8 × 10 −4 .We conclude that this minimal value is always very small compared to 1.We also remark that when Σ 2 0 goes to 1, i.e. for a virtually infinite dark photon mass m x , the ratio b/|a| is still very small since Finally, as Σ 2 0 decreases below Σ 2 min , the ratio b/|a| increases like ϵ 2 e 2 Q ′ 2 /12πΣ 2 0 and overcomes 1 below the critical value We are led to conclude that, as long as the mass degeneracy parameter Σ 2 0 is larger than Σ 2 cr , the ratio b/|a| is always less than 1.The thermal average ⟨σ ann v⟩ of the DM annihilation cross-section behaves like 1/Σ 3 for DM dispersion velocities Σ slightly above Σ 0 .Close to the peak, the integral J(a, b) may actually be approximated by as shown in Appendix A, and the integrand of J ann is proportional to This integrand is maximal for |a| = 3/2, i.e. for the DM dispersion velocity Σ M equal to 2/3 (m x /2m ϕ )Σ 0 , a value which we will approximate hereafter by Σ 0 since the dark photon mass m x is very close to 2m ϕ .To summarize, after freeze-out, DM cools down while Σ decreases.The cross-section ⟨σ ann v⟩ increases like 1/Σ 2 and, close to the peak where most of DM annihilation occurs, like 1/Σ 3 .At the lower boundary y M of the integral J ann , the DM dispersion velocity is Σ M ≃ Σ 0 and the DM temperature T ϕ is equal to m ϕ Σ 2 0 .
2. Case Σ 2 0 < Σ 2 cr -the non-resonant regime For a mass degeneracy parameter Σ 2 0 less than Σ 2 cr , the ratio b/|a| is this time larger than 1.The integral J(a, b) is equal to 1 as long as b is smaller than 1.As shown in Appendix A, beyond that point, J(a, b) can be approximated by (III.25) The transition occurs at b = √ 3/2, at a dispersion velocity Σ M which is equal to 2/ √ 3 1/2 (m x /2m ϕ ) Λ 0 .For simplicity, we will use from now on Λ 0 as the benchmark value for Σ M in this case.In the early Universe, after freeze-out, ⟨σ ann v⟩ increases like 1/Σ 2 until the DM dispersion velocity Σ has decreased down to Λ 0 .Below that point, the annihilation cross-section drops like Σ 2 and vanishes.The lower boundary y M of integral J ann is defined in such a way that the DM dispersion velocity is Σ M ≃ Λ 0 and the DM temperature T ϕ is equal to m ϕ Λ 2 0 .

C. Results without kinetic decoupling
In Fig. 3, we consider a benchmark example of a 1 GeV DM scalar particle.The dark charge g x and the mixing angle ϵ have been respectively set equal to 0.1 and 10 −6 .In the left panel, the green curve features the evolution of the DM relic abundance Ω ϕ h 2 as a function of the mass degeneracy parameter Σ 2 0 .This curve has been obtained numerically under the assumption that DM does not undergo kinetic decoupling.We have first determined the freezeout point x F by solving equation III.10 with the help of a dichotomy.We have then calculated I ann by integrating expression III.12.At this stage, we have assumed that DM has always the same temperature as the SM plasma.This assumption will be challenged in the forthcoming section III D. The relic codensity ñ0 has been derived from III.11, relations III.13 and III.14 yielding Ω ϕ h 2 .
A close inspection of the green curve allows us to clearly identify two distinct regimes.For values of Σ 2 0 smaller than approximately 10 −13 , Ω ϕ h 2 is constant and the curve exhibits a plateau.As Σ 2 0 increases above that value, the relic abundance first decreases, reaches a minimum slightly below 10 −6 and eventually increases to explode close to the boundary Σ 2 0 = 1.The curve exhibits a plateau for small values of Σ 2 0 and a trough for the larger values.Both features are characteristic and can be understood with the help of the approximations developed in Sec.III B. In particular, we will try to identify the values of Σ 2 cr and Σ 2 min along the green curve.As a general remark, we notice that since DM is assumed here to be always thermally coupled to the SM plasma, the temperatures of both components are equal at all times, hence the identity In the expression III.15 of J ann , the variable y can be identified with the DM dispersion velocity squared Σ 2 .This implies that the lower boundary y M is equal to Σ 2 M .

The plateau
For values of Σ 2 0 smaller than Σ 2 cr , the ratio b/|a| overcomes 1 (see Fig. 2).In this regime, the annihilation crosssection ⟨σ ann v⟩ is maximal for Σ M ≃ Λ 0 , and the integral J ann is performed from y M ≃ Λ 2 0 to y F .In this interval, the function J(a, b) is equal to 1 and we get Since the mass degeneracy parameter Σ 2 0 is smaller than Σ 2 cr , which is itself smaller than Σ 2 min , the reduced decay width simplifies to the constant Physically, the decay of a dark photon into a pair of DM scalars ϕ φ is kinematically suppressed when its mass is nearly degenerate with 2m ϕ .In these conditions, the dark photon decays only into fermion pairs.We notice that the reduced width Λ 2 0 , the integral J ann and eventually the relic abundance Ω ϕ h 2 no longer depend on Σ 2 0 , hence the plateau.The mass degeneracy parameter Σ 2 0 has disappeared from the problem.At this stage, we will not compare our approximation III.17 with the numerical value of Ω ϕ h 2 ≃ 2.5 × 10 4 read from Fig. 3, insofar as we are mostly concerned in this article with kinetic decoupling.There is also an additional complication.The plateau is expected to extend up to Σ 2 cr , which is equal to 9.7 × 10 −15 in our example.This value corresponds to a DM temperature, and therefore to a plasma temperature in the absence of kinetic decoupling, of 9.7 µeV.But today the CMB temperature T 0 is equal to 235 µeV.We are in the particular situation where y 0 is larger than y M ≃ Λ 2 0 ≡ Σ 2 cr .This has two consequences.The integral J ann must be performed from y 0 to y F , and not from y M to y F .This yields the result ln (T F /T 0 ) instead of III.27.Then, the plateau extends beyond the boundary Σ 2 cr , and this as long as the Breit-Wigner enhancement of the cross-section ⟨σ ann v⟩ does not perturb too much the integral J ann .At least, this is so as long as Σ 2 0 is less than y 0 , provided that the ratio b/|a| is not too small compared to 1, and that the function J(a, b) is close to unity.In our example, we anticipate an extension of the plateau up to y 0 = T 0 /m ϕ = 2.35 × 10 −13 , more than an order of magnitude above Σ 2 cr .In Fig. 3, the transition occurs around 4 × 10 −13 , not too far from what is expected.

The trough
For values of Σ 2 0 larger than Σ 2 cr , the ratio b/|a| is smaller than 1.The annihilation cross-section undergoes a Breit-Wigner enhancement which peaks at the dispersion velocity Σ M ≃ Σ 0 .The integral J ann is performed from y M ≃ Σ 2 0 to y F .Above the peak, ⟨σ ann v⟩ scales like 1/Σ 3 .If DM annihilation reaches its highest intensity before the present epoch, i.e. if y M is larger than y 0 , the integral J ann is significantly modified with respect to the case of the plateau studied above.Actually, for values of b/|a| sufficiently small with respect to 1, the function J(a, b) can be replaced, in the integrand of J ann , by its approximation III.23, hence The ratio |a|/b is by definition identical to Σ 2 0 /Λ 2 0 , while Σ 2 is equal to the parameter y insofar as DM is kinetically coupled to the SM plasma.We have dropped the exponential and replaced it by the sharp boundary at y M .Neglecting the term with y F , and noticing that m x and 2m ϕ are quasi-degenerate as long as Σ 2 0 is not too close to 1, we get Inserting this result into relation III.17 yields the DM relic abundance . (III.31) The mass degeneracy parameter Σ 2 0 appears in the argument T M ≃ m ϕ Σ 2 0 of the plasma function inside the factor F, and in the ratio Λ 2 0 /Σ 2 0 .In the former case, it has little impact on Ω ϕ h 2 since P is a slowly varying function of plasma temperature.The DM relic abundance depends on Σ 2 0 essentially through the ratio Λ 2 0 /Σ 2 0 .We can readily apply our analysis of Sec.III B. As Σ 2 0 increases from Σ 2 cr to the upper limit 1, the ratio Λ 2 0 /Σ 2 0 = b/|a|, and hence Ω ϕ h 2 to which it is proportional according to Eq. III.31, decreases to a minimum reached at Σ 2 min and then increases.This is actually what we observe in Fig. 3 where the green curve exhibits a trough.Its minimum corresponds to a DM dispersion velocity Σ 2 min of order 4×10 −7 .This value can be compared to our expectation III.19.Plugging in it the numerical values of g x and ϵ, and noticing that Q ′ 2 is very close to 4 for a 1 GeV DM candidate, we derive a value of 4.4 × 10 −7 , in excellent agreement with the numerical result.This gives us confidence in our approach and put our approximation on firm grounds.
According to Fig. 3, below Σ 2 min , the relic abundance Ω ϕ h 2 decreases with Σ 0 as a power law with index close to −2.Above the minimum, it increases like Σ 0 , following a power law with index +1.For Σ 2 0 larger than 10 −2 , Ω ϕ h 2 sharply increases.To understand the trough which the green curve exhibits, we can start from the approximation III.31 for Ω ϕ h 2 and from the expression III.18 of the ratio Λ 2 0 /Σ 2 0 .Below Σ 2 min , i.e. for small values of Σ 2 0 , that ratio simplifies to Noticing that the effective charges Q ′ 2 and Q2 are essentially equal as long as m x is nearly degenerate with 2m ϕ , we infer the relic abundance We do find that Ω ϕ h 2 scales as 1/Σ 2 0 , as observed in Fig. 3.We also notice that the relic abundance depends on the dark charge g x and no longer on the mixing angle ϵ.
For values of Σ 2 0 larger than Σ 2 min , the ratio Λ 2 0 /Σ 2 0 simplifies this time to g 2 x Σ 0 /48π.This yields Our analysis allows us to understand why the DM relic abundance follows a power law in Σ 0 with index +1 above Σ 2 min .We also remark that Ω ϕ h 2 depends in this case on the mixing angle ϵ and not on the dark charge g x .When Σ 2 0 is close to 1, our analysis and relation III.34 no longer apply.The dark photon mass m x starts to be much larger than 2m ϕ .Since the ratio b/|a| is still smaller than 1, a Breit-Wigner resonance may enhance the annihilation cross-section.But the peak dispersion velocity Σ M is larger than Σ 0 by a factor m x /2m ϕ which could be very large.If so, Σ M exceeds the speed of light and ⟨σ ann v⟩ is never enhanced, whatever the DM dispersion velocity.DM annihilation is, in this case, p-wave dominated.The dark photon propagator in the s-channel yields furthermore a factor 1/m 4 x which suppresses DM annihilation when m x diverges.It is no surprise then if Ω ϕ h 2 explodes when Σ 0 goes to 1, as observed in Fig. 3.
As a side remark, we notice that depending on the values of g x and ϵ, the dispersion velocity Σ 2 min may well exceed the boundary 1.This happens actually for small g x and large ϵ.In this case, the general behavior of Ω ϕ h 2 as a function of Σ 2 0 is not qualitatively changed with respect to what has been discussed above.We still have a plateau below Σ 2 cr and a trough above, the relic abundance starting to decrease, reaching a minimum and eventually sharply increasing close to 1.The position of the minimum is no longer given by Σ 2  min , but must be defined numerically.

D. Results in the presence of kinetic decoupling
We have so far assumed that DM is in thermal contact with the SM plasma, and that their respective temperatures T ϕ and T are equal at all times relevant to our analysis.We now challenge this assumption.Thermalization of DM occurs primarily through an exchange of energy due to collisions with the SM plasma.Should DM be slightly colder than the plasma, for instance, the latter would reheat the former by injecting energy through collisions of SM fermions on DM scalars.If this occurs fast enough, T ϕ relaxes rapidly toward T and DM is thermalized.In appendix B, we present a simplified treatment of that process.Starting from Lagrangian II.1, we derive the rate at which energy is exchanged between DM and the plasma.DM annihilations must also be included together with collisions to determine how both DM density and temperature vary concomitantly.In appendix B, we actually demonstrate that, under the specific assumptions presented at the beginning of this section, reaction III.1 is able alone to thermalize DM with the plasma, should collisions be turned off.Taking into account both annihilations and collisions, we establish the simplified differential equation B.27 which the DM temperature follows.This allows to define the rate Γ KD rel with which T ϕ relaxes toward T , and the rate Γ KD eq with which the kinetic equilibrium itself evolves.Kinetic decoupling occurs at temperature T KD , for which both rates are equal.
In the right panel of Fig. 3, the red curve features the evolution of Ω ϕ h 2 as a function of the mass degeneracy parameter Σ 2 0 when kinetic decoupling is taken into account.The other parameters have the same values as for the green curve.The red curve has been obtained numerically by first solving Eqs.III.10 and B.31 using dichotomies.This allows to determine the freeze-out T F and kinetic decoupling T KD temperatures.As explained in appendix B, kinetic decoupling occurs after freeze-out, hence the sequence T F ≥ T KD .We then calculate I ann by integrating expression III.12 from freeze-out until now.Between freeze-out and kinetic decoupling, i.e. as long as T is larger than With KD Without KD FIG.3: DM relic abundance as a function of Σ 2 0 for fixed gx and ϵ (respectively 10 −1 and 10 −6 ), and a DM mass of 1 GeV.In the left panel, the case without kinetic decoupling of the DM is highlighted in green with the asymptotic scalings around Σ 2 min .In the right panel, the case with kinetic decoupling is highlighted in red, with the asymptotic scalings around Σ2 min and Σ 2 cr .The horizontal black line indicates the Planck DM relic abundance ΩDMh 2 .
T KD , both DM and plasma temperatures are equal.After kinetic decoupling, DM behaves as a non-relativistic gas undergoing adiabatic cooling and its temperature T ϕ drops as a −2 , with a the scale factor of the Universe.The plasma also cools down adiabatically, but is now decoupled from DM.The temperatures of both components are related by We then determine the relic codensity ñ0 from III.11 and derive the DM relic abundance Ω ϕ h 2 using relations III.13 and III.14.The red and green curves of Fig. 3 are qualitatively similar.Both exhibit a plateau for small values of Σ 2 0 and a trough for larger values.The red curve becomes flat below 5 × 10 −15 .For larger values of Σ 2 0 , the DM relic abundance starts to decrease and reaches a minimum for a mass degeneracy parameter of order 2 × 10 −7 , slightly below Σ 2 min .The red curve then follows a rising power law and eventually sharply increases close to 1.The essential effect of including kinetic decoupling in the calculation of Ω ϕ h 2 is to shift the curve downward and to get smaller values of the DM relic abundance.This can be easily understood.When DM decouples thermally from the primordial plasma, its temperature drops faster than if thermal contact was continuously established.We have just showed that T ϕ decreases as a −2 , while T scales approximately like a −1 , where a is the scale factor of the Universe.As DM cools down, the annihilation cross-section ⟨σ ann v⟩ increases.It peaks at the DM dispersion velocity Σ M , where most of the annihilation takes place.When kinetic decoupling is included, this occurs at a higher plasma temperature T M , i.e. at an earlier time t M when the DM population is denser.A stronger DM annihilation at ⟨σ ann v⟩ peak results in a smaller relic abundance, hence the observed shift between the green and the red curves.When kinetic decoupling is included, the relation between Σ 2 and y is also slightly more involved.In Sec.III C, we could identify y and Σ 2 at all times.This is only possible now before kinetic decoupling.For temperatures below T KD , the new relation is where y KD denotes the ratio T KD /m ϕ .If most of DM annihilation occurs after kinetic decoupling, i.e. if T M is well below T KD , the integral J ann may be approximated by

37)
The ratio h eff (T KD )/h eff (T ) varies slowly in time and we have taken it at peak annihilation when the plasma temperature is T M .We also remark that the relation between y M and the DM dispersion velocity Σ M has become Equipped with these notations, we are ready to analyze the red curve.

The plateau
As showed in Sec.III C, the ratio b/|a| overcomes 1 when Σ 2 0 is smaller than Σ 2 cr .This implies a transition at 9.7 × 10 −15 , in agreement with the value of 5 × 10 −15 mentioned above.In the plateau regime, the annihilation cross-section ⟨σ ann v⟩ is maximal for Σ M ≃ Λ 0 and we get Using this expression into relation III.17 leads to the approximate DM relic abundance where F KD denotes the pre-factor Along the plateau, i.e. below Σ 2 cr , the dark photon mass m x is very close to 2m ϕ and we can identify Q with Q ′ .We have also denoted by x KD the inverse of y KD .We notice that the mass degeneracy parameter Σ 2 0 has disappeared from the problem.With the values of model parameters considered in Fig. 3, we find that along the plateau, freeze-out and kinetic decoupling occur respectively at x F = 1.71 and x KD = 8.48.This translates into a decoupling temperature T KD of 0.118 GeV.The peak of DM annihilation is reached at the plasma temperature T M = 58 eV well above T 0 .Relation III.40 yields a DM relic abundance of 0.139, to be compared to the numerical value of 0.220, obtained asymptotically by setting the mass degeneracy parameter Σ 2 0 at 10 −17 .This result is very encouraging given the simplicity of our approximation.We could improve the agreement by slightly increasing the critical DM dispersion velocity Σ M above Λ 0 .Shifting it from Λ 0 to 3Λ 0 /2, for instance, yields a relic abundance of 0.208, in better agreement with the numerical result.

The trough
For values of Σ 2 0 larger than Σ 2 cr , the ratio b/|a| is smaller than 1 and the analysis proceeds along the same line as in Sec.III C. Thermal decoupling complicates the relation between the y parameter and the DM dispersion velocity Σ.In the trough regime, DM annihilation is Breit-Wigner enhanced at its peak.Replacing in the integrand of expression III.37 the function J(a, b) by its approximation III.23 leads, after some algebra, to the integral Making use once again of the conversion III.36 and setting the peak dispersion velocity Σ M at Σ 0 , we get As long as Σ 2 0 is not too close to 1, we can identify m x with 2m ϕ .Combining relations III.43 and III.17 leads to the DM relic abundance .

(III.44)
To test this approximation, we have varied Σ 2 0 from 10 −13 up to 0.1, and compared the numerical result with the value given by III.44.On a vast region of that interval, both results agree at the percent level and, in some cases, even at the per mille level.The agreement lessens close to the minimum located at Σ 2 0 ≃ 1.66 × 10 −7 , where the approximation yields a relic abundance of 6.12 × 10 −5 to be compared to the numerical result of 7.45 × 10 −5 .We finally notice that above Σ 2 0 ≃ 5 × 10 −3 , the numerical relic abundance increases more sharply than its approximation III.44.Reasons for this have already been presented in Sec.III C. Along the trough, Ω ϕ h 2 is proportional to the ratio which is minimum at the mass degeneracy parameter So is Ω ϕ h 2 .If thermal contact was constantly established between DM and the SM plasma, the DM relic abundance would exhibit a minimum at Σ 2 min .In the presence of kinetic decoupling, the new minimum is a factor 4 2/3 smaller.In Fig. 3, this shift of the minimum of Ω ϕ h 2 between the green and red curves can be clearly seen.Relation III.46 yields a value of 1.75 × 10 −7 for the minimum of the red curve, close to the numerical value of 1.66 × 10 −7 .
Below Σ2 min , the ratio Λ 2 0 /Σ 0 can be approximated by ϵ 2 e 2 Q ′ 2 /12πΣ 0 and expression III.44 simplifies into where we have identified Q2 with Q ′ 2 .Above the minimum, we can replace Λ 2 0 /Σ 0 by (g 2 x /48π)Σ 2 0 to get In Fig. 3, the red curve decreases actually with Σ 0 as a power law of index −1 while, above the minimum, it follows a power law with index +2.We also remark that Ω ϕ h 2 depends in the former case on the dark charge g x and not on the mixing angle ϵ, while it is the opposite in the latter case.This property will play a crucial role in Sec.IV and will help understand the results.For completeness, we have derived an approximate expression for the minimal DM relic abundance.Since the minimal value of the ratio Λ 2 0 /Σ 0 is given by we infer that the DM relic abundance, in the presence of kinetic decoupling, reaches a minimal value of where we have identified once again Q2 with Q ′ 2 .This expression yields a minimum value of 6.13×10 −5 to be compared to the numerical result of 7.46 × 10 −5 for a mass degeneracy parameter Σ 2 0 of 1.75 × 10 −7 .Since approximation III.50 is based only on the behavior of the ratio Λ 2 0 /Σ 0 , it disregards the influence of the temperatures T KD and T M on Ω ϕ h 2 , hence the very small differences with the results quoted above.For Σ 2 0 larger than 10 −2 , Ω ϕ h 2 sharply increases.In the next section, we will impose that the DM relic abundance is equal to the observed value Ω DM h 2 of 0.1200 [1].For this, at fixed DM mass m ϕ , dark charge g x and mixing angle ϵ, we will look for values of the mass degeneracy parameter Σ 2 0 that fulfill this requirement.Depending on the height of the red curve with respect to the level of Ω DM h 2 , three configurations are possible.
1.If the red curve is too high, there is no solution and the DM relic abundance Ω ϕ h 2 is always larger than Ω DM h 2 .
The Universe is overclosed by scalar DM.
2. In the configurations of interest, two values of Σ 2 0 fulfill the condition Ω ϕ h 2 = Ω DM h 2 .One solution lies on the decreasing left branch of the red curve, below the critical value Σ2 min at which the DM relic abundance is minimal.The other solution is located above Σ2 min , on the rising right branch of the red curve.These solutions are dubbed hereafter left and right branch solutions.
3. If the red curve is too low, there is no left branch solution insofar as the plateau stands below Ω DM h 2 .The Universe is underclosed by scalar DM below Σ2 min .A right branch solution still exists.

IV. RESULTS
We present our results in the form of scans on the DM model parameter space in the g x and ϵ plane for three chosen DM particle masses: 200 MeV, 1 GeV and 5 GeV.We take into account kinetic decoupling as presented in Sect.III D. At 200 MeV, DM annihilation is dominated by the leptonic channels whereas at 1 and 5 GeV, there is also a significant contribution from DM annihilating directly into quark pairs.Moreover, for these masses, we mostly avoid the stringent direct detection limits.The impact of various constraints on the parameter space is summarized in Figs. 6, 7, and 8, both for the left and right branch cases.In Sect.IV A, we first discuss how the DM relic abundance shapes the allowed region in the g x and ϵ plane.Then, in Sect.IV B, IV C, and IV D, we show how other constraints, from DM annihilation in the Milky Way, in the early Universe, and from direct detection and accelerators, shrink the allowed region of the DM parameter space.
A. Limits on gx and ϵ set by the relic abundance ΩDMh 2 As discussed previously in Sect.III D notably when commenting Fig. 3, requiring Ω ϕ h 2 = Ω DM h 2 leads, at fixed values of the parameters g x and ϵ, to two different solutions, depending on the mass degeneracy parameter Σ 2 0 .We first discuss the solutions where Σ 2 0 < Σ2 min which we dubbed left branch solutions in the previous section, and then consider the right branch solutions corresponding to Σ 2 0 > Σ2 min .
1. Left branch solutions -Σ 2 0 < Σ2 min Figure 4 shows the allowed regions for the case of a 1 GeV DM scalar after enforcing Ω ϕ h 2 = Ω DM h 2 with various quantities displayed in the color bar.The trends explained in the following are similar for the other DM masses.The top panels of Fig. 4 feature in their respective color bars the kinetic decoupling parameter x KD (left) and the mass degeneracy parameter Σ 2 0 (right).In the colored band, the relic density constraint is satisfied, while in the white regions above and below either Ω ϕ h 2 < Ω DM h 2 or Ω ϕ h 2 > Ω DM h 2 .First, in the region where Ω ϕ h 2 < Ω DM h 2 , the value of the plateau of Ω ϕ h 2 at low Σ 2 0 (see right panel of Fig. 3) does not exceed Ω DM h 2 .Since the plateau height scales as g −2 x ϵ −1 , at the boundary where Ω ϕ h 2 = Ω DM h 2 , we obtain the following scaling: ϵ ∝ g −2 x , which is the slope observed in the upper-right corner of the figures, between the white and colored regions.Second, in the region where Ω ϕ h 2 > Ω DM h 2 , the minimal value of Ω ϕ h 2 at Σ 2 0 = Σ2 min is always larger than Ω DM h 2 .Given the scaling of Ω ϕ h 2 min shown in Eq.III.50, at the boundary where Ω ϕ h 2 min = Ω DM h 2 , we recover the scaling ϵ ∝ g −2 x .This slope is observed in the lower-left corner, between the white and colored regions in the figures.
Note that the discontinuity observed around g x = 10 −3 stems from the QCD phase transition as we shall explain hereafter.As shown in appendix B 1, kinetic decoupling happens after freeze-out, and even possibly after the QCD phase transition.In the black region of the upper-left panel of Fig. 4, x KD ≈ 1, which means that, since m ϕ = 1 GeV, T KD ≈ 1 GeV.This is a temperature above the QCD phase transition.In the green region, x KD ≈ 10, thus T KD ≈ 0.1 GeV which is a temperature below the QCD phase transition.As explained in Sec.III D, the later the kinetic decoupling, the larger the relic abundance.Thus, for given model parameters, if T KD becomes smaller than T QCD , the temperature T ϕ starts dropping later, in particular, because of latent heat release, and annihilation occurs in a less dense medium resulting in a larger value for Ω ϕ h 2 min .Because the latter scales as ϵ −2/3 , one has to boost the value of ϵ to satisfy Ω ϕ h 2 = Ω DM h 2 , hence the jump in the plot observed around g x ≈ 10 −3 .
Another feature we would like to explain is the sharp vertical cut in the allowed parameter space around g x ≈ 10 −4 .In this region of the (g x , ϵ) plane, the value of Σ2 min , as given by Eq.III.46, becomes larger than 1.In the trough regime, as the mass degeneracy parameter Σ 2 0 increases, the relic abundance decreases according to relation III.47.
A minimum is reached near Σ 2 0 ≃ 0.01.This corresponds here to a plasma temperature T M ≃ 10 MeV.Requiring that this minimum is less than Ω DM h 2 sets a lower limit on g x .With a kinetic decoupling temperature T KD of order 50 MeV, using relation III.47 leads to an approximate value for that bound of g x = 1.4 × 10 −4 , in excellent agreement with the numerical result.We reproduced this exercise for a DM mass of 200 MeV and 5 GeV, and found limiting g x values of 3.0 × 10 −5 and 5.2 × 10 −4 , in agreement with the summary plots Fig. 6 and Fig. 8, respectively.The remaining panels in Fig. 4 will be discussed in the following subsections.The corresponding figures for the right branch solutions are shown in Fig. 5.The shape of the lower boundary, below which Ω ϕ h 2 > Ω DM h 2 is exactly the same as in the left branch case, and is explained by the same arguments.The main difference with the left branch solution is that the upper boundary above which Ω ϕ h 2 < Ω DM h 2 does not exist in this case.In the right branch, it is always possible to increase the value of Σ 2 0 and thus to increase Ω ϕ h 2 until we reach the required value.When Σ 2 0 approaches 1, the dark photon mass m x explodes and the Breit-Wigner resonance disappears.

B. Limits on gx and ϵ from the annihilation cross-section today
As discussed in Sect.II, the peculiarity of our model is the non-trivial dependence of the dark matter annihilation cross-section with the dispersion velocity Σ 2 , which, in the Breit-Wigner regime, peaks at Σ 2 0 .The value of Σ 2 0 is set by the requirement Ω ϕ h 2 = Ω DM h 2 .Using the scalings of Eq.III.47 and Eq.III.48 (also shown in the right panel of Fig. 3) implies that Σ 2 0 scales as g −4 x for the left branch and as ϵ 2 for the right branch.The resulting Σ 2 0 values are shown with a color scale, for these two cases, in the top-right panels of Fig. 4 and Fig. 5.As predicted by the scalings, we note the almost complete independence of Σ 2 0 values on ϵ and g x in the first and second case, respectively.From the values of Σ 2 0 , g x and ϵ, we compute the corresponding DM annihilation cross-section today in the Milky Way ⟨σ ann v⟩ MW in order to compare with observational limits and to set constraints on our model.For this, we use relation II.10 and set the DM dispersion velocity Σ equal to its value Σ MW in the Milky Way halo.Strictly speaking, the chosen value for Σ MW should depend on the position in the Milky Way (see e.g.[41]).However, given that the X-ray constraints that we will be using are drawn from multiple lines of sight, for simplicity, we set Σ 2 MW to the fiducial value 3 × 10 −7 c 2 .This value is obtained for an isothermal DM halo accounting for a flat rotation curve with the typical 230 km/s circular velocity.We anticipate that ⟨σ ann v⟩ MW is maximal at the Breit-Wigner peak, when the Milky Way dispersion velocity Σ 2 MW is of order the mass degeneracy parameter Σ 2 0 .We display ⟨σ ann v⟩ MW with a color scale in the lower-left panels of Fig. 4 and Fig. 5, for the left and right branches, respectively.Hereafter, we discuss these two cases independently.
In the lower-left panel of Fig. 4, the color shows that the DM annihilation cross-section values ⟨σ ann v⟩ MW are almost independent of ϵ.Going from small to large g x values, ⟨σ ann v⟩ MW increases rapidly, reaching a maximum at the Breit-Wigner peak, followed by a slow decrease toward the largest g x values.In fact, this behavior can be easily explained by the variations of ⟨σ ann v⟩ with Σ 2 as shown in Fig. 1.Going from small to large g x values implies decreasing Σ 2 0 values as shown in the upper-right panel.When g x increases, the Breit-Wigner peak is shifted toward lower DM dispersion velocities.The annihilation cross-section in the Milky Way follows a p-wave behavior and increases until a maximum is reached when Σ 2 MW is of order Σ 2 0 .This occurs at g x of order 10 −3 , above which ⟨σ ann v⟩ MW decreases.The reasoning developed above is only valid at first order, since we notice deviations from a unique dependence of ⟨σ ann v⟩ MW on g x , at large enough ϵ values.For example, we see that for a given value of g x , e.g. 10 −2 , going from small to large ϵ values, ⟨σ ann v⟩ MW is constant and then increases.This leads to a chevron-like feature for ⟨σ ann v⟩ MW in the (g x , ϵ) plane.We explain this feature by the change of behavior of ⟨σ ann v⟩ going from the Breit-Wigner resonance to the high-velocity regime.From Eq. II.10, we know that this transition typically occurs when J(a, b) starts saturating to 1 for low values of |a| (see also Fig. 9).When this happens, ⟨σ ann v⟩ MW ∝ g 2 x ϵ 2 , a scaling which is observed in Fig. 4. When the parameter |a| is not too small, the function J(a, b) is given by its approximation A.3.The DM annihilation cross-section is proportional to Σ 2 0 /Λ 2 0 .Along the left branch, we can use relation III.32 and we find that ⟨σ ann v⟩ scales like g 2 x Σ 3 0 .The tip of the chevron, i.e. the position of this transition, can be found by equating the approximation J 1 (a, b), valid at the Breit-Wigner peak, to 1. From Eq. A.3, this means that √ πa(|a|/b) = 1, and we readily get that ϵ 2 ∝ Σ 3 0 .Hence, as Σ 2 0 ∝ g −4 x for left branch solutions, we find that ϵ ∝ Σ 3/2 0 ∝ g −3 x , which corresponds to the line one can draw through the tips of the chevrons in the lower-left panel.

Right branch solutions -
On the lower-left panel of Fig. 5, we see that ⟨σ ann v⟩ MW values are almost independent of g x .Going from large to small ϵ values, this cross-section exhibits a rapid increase, reaches a maximum, and then decreases with ϵ.In the same way, as for the left branch case, this behavior is explained by the variations of ⟨σ ann v⟩ with Σ 2 as shown in Fig. 1.Going from large to small ϵ values implies decreasing Σ 2 0 values (upper-right panel of Fig. 5), and shifting the position of the Breit-Wigner peak of ⟨σ ann v⟩ toward smaller Σ 2 .Since ⟨σ ann v⟩ MW is evaluated at fixed DM dispersion velocity Σ 2 MW , the annihilation cross-section in the Milky Way follows a p-wave behavior as ϵ decreases.It increases until a maximum is reached when Σ 2 MW is of order Σ 2 0 .This occurs at ϵ of order 10 −8 .Below that value, ⟨σ ann v⟩ MW decreases.
As for the left branch, this picture is only valid at first order, important corrections occurring for ϵ > 10 −7 .For example, at g x = 10 −2 , going from ϵ = 10 −7 to 10 −2 , ⟨σ ann v⟩ MW decreases and then increases again.The latter behavior can be explained by going into some more detail.From Eq. II.10, we know that this behavior occurs in the p-wave regime when J(a, b) is well approximated by J 2 (a, b) (see Eq. A.1 and Fig. 9 of the appendix).In this approximation, one can easily show that ⟨σ ann v⟩ MW ∝ ϵ 2 /Σ 4 0 and, given the scaling of the right branch Σ 2 0 ∝ ϵ 2 (see Eq. III.48), we recover that ⟨σ ann v⟩ MW ∝ ϵ −2 .However, as one can notice from the right panel of Fig. 3, the scaling drawn from Eq. III.48 fails to reproduce the sharp rise of Ω ϕ h 2 as Σ 2 0 tends to 1. Instead, we can assume that Σ 2β 0 ∝ ϵ 2 with β being a number larger than 1.In that case ⟨σ ann v⟩ MW ∝ ϵ 2−4/β , which means that, when β gets larger than 2, ⟨σ ann v⟩ MW starts increasing with growing ϵ values.This explains the violet-blue spot of the lower-left panel of Fig. 5.

X-rays constraints
Light stable fermions (typically e + e − or µ + µ − ) from DM annihilating today, may energize the low-energy photons of the interstellar radiation field of the Galaxy and generate a sizable X-ray emission, via the Inverse Compton effect.In the lower-left panels of Fig. 4 and Fig. 5, we draw with a hatched region, the X-ray constraints of XMM Newton taken from [20], that we find to be the strongest of the literature.The constraints shown include only DM annihilating into e + e − pairs.This annihilation channel provides the strongest limit on ⟨σ ann v⟩ MW .The analysis by [20] provides the following upper bound on ⟨σ ann v⟩ MW for DM annihilating into e + e − : 7.8 × 10 −29 , 1.5 × 10 −28 , and 2.3 × 10 −27 cm 3 s −1 for 200 MeV, 1 GeV and 5 GeV DM mass, respectively.In our case, we apply these limits on the product B e + e − × ⟨σ ann v⟩ MW , with the branching ratio B e + e − ≡ Q2 e + e − / Q2 (see Eq. II.6, where Q2 e + e − only includes e + e − in the sum).The other constraints discussed in [20], which are based on µ + µ − and π + π − , are not competitive and have not been implemented.
C. Limits on gx and ϵ from CMB spectral and angular distortions DM annihilation in the early Universe injects energy in the primordial plasma and may generate distortions in the radiation spectrum.These distortions get frozen after recombination and leave indelible deviations of the CMB from a pure black body spectrum (see e.g.[24] for an introduction).The characteristics of the distortions can be addressed very generally by solving the so-called Kompaneets equation [42] which describes the Comptonization of photons by free thermal electrons.However, depending on the epoch at which the energy is released, approximations can be accurately used and avoid solving this non-linear equation.If the injection occurs for redshifts larger than z DC = 1.98 × 10 6 [43], the energy is rapidly redistributed to the photons via Compton scattering, while the number of photons is adjusted by photon non-conserving processes: double Compton (DC), which sets z DC , and thermal Bremsstrahlung.In this case, the CMB spectrum just undergoes a temperature shift.If the energy injection occurs after z DC , but before z C = 5.8 × 10 4 [43], the number of photons remains unchanged while Comptonization, which sets z C , ensures an efficient redistribution of energy among photons.The energy per photon is increased and leads to a non-vanishing chemical potential µ in the Bose-Einstein spectrum, referred to as µ-distortion.Finally, if the energy injection happens for redshifts smaller than z C , a typical Compton y-distortion arises when scatterings become inefficient in exchanging energy.This classification is useful to easily compare with present constraints from FIRAS [25,26] which set the upper limits |µ| < 4.7 × 10 −5 and |y| < 1.5 × 10 −5 , but in general, there is a large variety of possible distortions with a smooth transition from y to µ distortions that could give complementary information.In practice we compute the µ and y-distortions using the prescription presented in [43] and [44], such as: FIG. 5: Allowed parameter space in the (gx, ϵ) plane, for right branch solutions (Σ 2 0 > Σ2 min ) and for a DM mass of 1 GeV.From the upper-left to the lower-right panel, the color code shows xKD, Σ 2 0 , the corresponding DM annihilation cross-section in the Milky Way ⟨σannv⟩MW and the CMB µ-distortion, respectively.In the bottom panels, the hatched regions draw the exclusion constraints by XMM X-ray measurements and FIRAS (COBE) on the µ-distortion (solid line) and on the y-distortion (dashed line).The full interpretation of these plots is provided in the text.
where we have introduced the window functions J µ (z) and J y (z) defined by: (IV.4) We compute the energy released in the photon bath dρ γ = f em × ⟨σ ann v⟩n 2 ϕ dt × 2m ϕ corresponding to the energy injected by DM annihilation during dt.The pre-factor f em accounts for the fraction of the energy injected in the form of electromagnetic energy and depends on the dark matter mass.This factor has already been computed for this model in Ref. [30] (see their Eq.3.4 and Fig. 2), and we extract the following values: 0.52, 0.28, 0.27, for m ϕ = 200 MeV, 1 GeV and 5 GeV.However, in this paper, the authors were interested in the impact of the ionization on the CMB and so included ionization efficiencies.Hence, using their values for f em , we slightly underestimate the µ and y values we compute.The constraints we will draw from FIRAS [25,26] upper bounds on µ an y, thus give conservative limits in the (g x , ϵ) parameter space.
The results for the µ-distortion are shown with a color scale on the lower-right panel of Fig. 4 and Fig. 5, for the left and right branch, respectively.As for ⟨σ ann v⟩ MW , the µ values are primarily dependent on g x (resp.ϵ) in the left (resp.right) branch case.The evolution of µ with these parameters follows the same trend as for ⟨σ ann v⟩ MW .This is not a coincidence since, by definition, µ traces the DM annihilation history in the redshift window [z DC , z C ].In particular, we remark first that the peak of µ values occurs at larger g x (resp.smaller ϵ) compared to the peak in ⟨σ ann v⟩ MW in the left (resp.right) branch.This is explained by the fact that the DM dispersion velocity Σ 2 at the epoch [z DC , z C ] is smaller than the virialized one in the DM galactic halo today.Second, the peak of the µ values is broader (the red-colored region is wider) than the one seen on ⟨σ ann v⟩ MW .This is because the DM annihilation is integrated over the full redshift window [z DC , z C ], including a broad range of Σ 2 values with respect to the plot of ⟨σ ann v⟩ MW which displays a picture of the annihilation for the specific Σ 2 MW value.The hatched-region delineated by the solid line corresponds to the region of the parameter space excluded by the upper bound |µ| < 4.7 × 10 −5 from FIRAS [25,26].The results for the y-distortion look very similar in shape to the ones for the µ-distortion, and we only show with a dashed line the region of the parameter space excluded by the upper bound |y| < 1.5 × 10 −5 from FIRAS [25].We note that the constraints from the y-distortion are always weaker than the ones from the µ-distortion.In the summary plots of Fig. 6, 7 and 8, we call CMB constraints the limits from µ-distortion.
For completeness, we also report the constraints coming from CMB anisotropies.Following Ref. [30], we have computed the effective annihilation parameter defined in their Eq.(3.1): where R ≡ n ϕ (CMB)/n ϕ (today) being the ratio between the dark matter density at the CMB epoch and today.As the injection of energy in the primordial plasma from DM annihilation has maximal effect for z ≈ 600, we have computed p ann at that redshift.We have applied the limits on p ann from [1].Thus, for the left branch solutions (Σ 2 0 < Σ2 min ), the corresponding exclusion region stands to the right of the dotted line in the lower-right panel of Fig. 4 for a 1 GeV DM scalar.We find the same trend for the 200 MeV and 5 GeV cases, with much weaker constraints for the latter case.This is expected insofar as the peak of DM annihilation takes place earlier for heavier masses.For all the masses considered, we always find the p ann constraint to be sub-leading with respect to the µ-distortions one.For right branch solutions (Σ 2 0 > Σ2 min ), we find no constraints from p ann since the peak of DM annihilation occurs well before recombination.

D. Limits on gx and ϵ from direct detection and collider experiments
Other significant limits on the (g x , ϵ) parameter space come from the DM direct detection experiments as well as from the accelerator searches.

Direct detection limits
Dark photon interacting with ϕ and the SM particles can be constrained by the direct detection limits because DM can scatter off nucleons or electrons in the detector material through A ′µ exchange in ϕ SM → ϕ SM processes.The spin-independent elastic scattering cross-section is given by [45,46] where µ T is the reduced mass of DM with the target species T = electron/nucleon and q T is the U ′ (1) charge of the target particle, which is 1/2 for a proton, −1/2 for a neutron and 1 for an electron.q is the momentum transfer.For our parameter range of interest, m x ≈ 2m ϕ , which implies the zero-momentum transfer limit for the case of a heavy mediator, i.e. m x ≫ q.Eq.IV.6 can be approximated as: The limits from electron recoil experiments are relevant for lighter DM whereas for larger m ϕ ( ie, m ϕ ≳ 100 MeV), the limits from DM-nucleon scattering kick in.Although detection sensitivity rapidly decreases at low recoil energies, there have recently been substantial efforts in low-mass detection techniques, with experiments like PandaX-4T and XENON IT using liquid Xenon (LXe) detectors.DM scattering off the target nuclei produces scintillation photons (S1) and ionized electrons (S2) (produced through Migdal effect).Both S1 and S2 signals are used for low-mass DM detection, although the S2-only signal is more efficient toward lower energy.For nuclear recoil experiments, around ∼ 100 MeV DM mass, the most efficient limits come from the S2-only bounds.In this category, PandaX-4T limits improve the previous limits from CRESST-III and XENON-IT(M) by several orders [47,48].For heavier DM of a few GeV mass, the most stringent limits come from PandaX-4T (S1-S2), which improves the old XENON-IT(S2) limits by ∼ a factor of 2 [47,48].In the case of electron recoil, for DM mass m ϕ ≲ 10 MeV, the most stringent limits come from SENSEI [49], while for larger DM mass these limits are surpassed by DarkSide-50 [50] and eventually XENON-IT(S2) limits take over around m ϕ ∼ 30 MeV [51].Recent results from PandaX-4T improve XENON-1T limits by almost an order for DM mass ≳ 50 MeV [52].
In Figs.6-8, the exclusion limits by the existing direct detection experiments are illustrated in the (g x , ϵ) plane.As obvious from Eq. IV.7, the scattering cross-sections are effectively independent of Σ 2 0 (it only introduces a negligible correction when taking 2m ϕ ≃ m x ), which implies that the exclusion regions in (g x , ϵ) plane do not change for the left and right branch solutions.For m ϕ = 200 MeV, the orange shaded region in Fig. 6 corresponds to the region excluded by PandaX-4T (S2 only+Migdal) limits, which is the most constraining so far with the allowed upper limit of spinindependent scattering cross-section reaching down to σ SI N ≃ 3.25 × 10 −38 cm 2 [47].The electron recoil limits can be obtained from PandaX-4T (constant W model) [52], where the exclusion region corresponds to σ SI e ≳ 2.1 × 10 −41 cm 2 .These limits are much less constraining in comparison with the nuclear recoil bounds and therefore are not shown in the figure.In Fig 7 and 8, the orange colored exclusion regions are obtained from PandaX-4T (S2 only+Migdal) and PandaX-4T (S1-S2) [47] respectively.We quote the following numbers for the respective allowed upper limits on σ SI N : σ SI N ≃ 1.61 × 10 −39 cm 2 for m ϕ = 1 GeV and σ SI N ≃ 1.22 × 10 −44 cm 2 for m ϕ = 5 GeV.We find that the electron recoil limits are practically irrelevant for these masses and do not constrain the range of g x and ϵ considered in the figures.
In Fig. 6 -8, the projected limits from the near-future direct detection experiments are shown with orange dashed lines.For 200 MeV and 1 GeV DM, we show projections from DARKSPHERE [53], an experiment proposed by the NEWS-G collaboration.DARKSPHERE uses a spherical proportional counter, optimised for detecting nuclear recoils with sub-keV energy.We quote the projected exclusions of σ SI N ≳ 2.47 × 10 −42 cm 2 for m ϕ = 200 MeV and σ SI N ≳ 1.7 × 10 −43 cm 2 for m ϕ = 1 GeV.For 5 GeV DM, the best projections are obtained from the SBC experiment [54,55] which uses liquid Argon (LAr) spiked with liquid Xenon (LXe) in a bubble chamber setup.We have used the projected sensitivity of σ SI N ≃ 7.28 × 10 −46 cm 2 for m ϕ = 5 GeV in the figures.

Accelerator limits
Dark photons can also be probed in various accelerator searches.Depending on the mass scale, the production and the decay of the dark photon in these experiments dictate the sensitivity of detection.While sub-GeV dark photons are best probed in the electron and the proton beam-dumps and other fixed-target experiments, they can also be tested in e + e − colliders, up to a few GeV mass [36,56].For dark photons of mass ≳ 10 GeV, there are constraints from several LHC searches [57,58].The signal for all these probes can come from both visible and invisible decay of dark photons.For visible dark photon searches, the typical signal is a pair of leptons whereas the invisible searches are sensitive to the missing energy signal, due to the dominant decay of the dark photon into DM.
In our model, the bounds can be obtained from both visible and invisible dark photon accelerator searches.The partial decay widths of the dark photon into a pair of DM and SM particles respectively can be obtained from Eq. II.14, which implies that BR(A When both decay modes are kinematically accessible, one should take into account the respective branching ratios while computing the constraints.The constraints will therefore depend on g x , Σ 0 and ϵ.The constraints on dark photons found in the literature only depend on ϵ since it is generally assumed that the branching ratio for either visible or invisible decay modes is 1.The ϵ dependence comes only from the dark photon production which is determined by its interactions with SM particles.For the dark photon mass range of our interest, the following search strategies yield significant limits: 1. visible searches : The best limits for dark photons of mass ∼ 100 MeV are obtained from BaBar [59] in e + e − collisions.The production channel for the dark photon is e + e − → γ A ′ , and the signal is observed when A ′µ promptly decays into visible final states: A ′ → e + e − and µ + µ − .BaBar provides the most stringent limits on visible dark photon decays for 100 MeV ≲ m x ≲ 200 MeV and 1 GeV ≲ m x ≲ 10 GeV.For the mass window 2 m µ ≲ m x ≲ 0.5 GeV, the prompt and the displaced vertex searches (A ′ → µ + µ − ) at LHCb produce the best limits [57].Here A ′µ is produced through meson decays : GeV is best constrained by the di-muon searches at LHCb and CMS [60].A heavier dark photon, which is relevant for these searches can be produced in Drell-Yann process qq ′ → A ′ at the LHC.
2. invisible searches : Light dark photons decaying into a pair of DM particles can be probed in the invisible search experiments.For m x in our range of interest, there are only a few experiments that produce the relevant constraints.For both Babar and LEP, the process under consideration is the production of an ordinary photon accompanied by a dark photon.The decay channel A ′ → ϕϕ † gives a signature of monophoton and missing energy.For 25 MeV ≲ m x ≲ 8 GeV, BaBar limits are the most stringent [61].For m x ≳ 8 GeV, LEP limits prevail [56,62,63].
In Figs.6-8, the exclusion regions corresponding to the accelerator limits are shown in magenta shade.We find that the dark photon decays almost entirely into SM particles in the parameter space of interest for the left branch solutions, therefore we use only the limits from visible decays.For the right branch solutions, the dark photon can either decay into SM particles or mostly invisibly, the latter occurs typically in the region at large g x .Therefore, we take the limits from both visible and invisible decay and rescale them in proportion to their respective branching ratios.For m ϕ = 200 MeV, (which implies m x = 400 MeV), the LHCb di-muon searches constrain ϵ LHCb ∼ 5 × 10 −4 .In addition, the displaced vertex searches from the proton beam-dump such as ν − CAL [36,56,64] and CHARM [65] constrain smaller values of ϵ, namely, 10 −7 ≲ ϵ ≲ 9 × 10 −7 .In these searches, the dark photon is produced from the meson decays such as η (η ′ ) → γ A ′ and is subsequently decayed within the detector into a pair of displaced muons.As mentioned above, the exclusion bands in Fig. 6 (left) correspond to the visible decay limits alone.In Fig. 6 (right), the magenta-shaded upper band corresponds to the region excluded by LHCb (for visible decay) and BaBar (for invisible decays).When the branching ratios are ∼ 100%, the limits on ϵ correspond to ϵ LHCb ∼ 5 × 10 −4 and ϵ BaBar ∼ 1 × 10 −3 .As obvious from the expressions for the branching ratios above, we find that toward smaller g x and larger ϵ, the visible decays dominate and as g x increases, invisible decays gradually take over.Thus the magenta band for smaller g x corresponds to the purely visible decay limits from LHCb and large g x to the purely invisible decay bounds from BaBar.For intermediate g x , where we find substantial contribution for both decay modes, the limits are rescaled with the respective branching ratios as the number of signal events is ∝ ϵ 2 × BR.However, for the lower band, the only relevant limits are for visible decays.Therefore, for the right branch solution, the bound gradually weakens for larger g x where we find a dominant invisible decay contribution.Fig. 7 corresponds to m x ∼ 2 GeV, which is best constrained by BaBar for both visible and invisible searches.The allowed upper limits of ϵ from both types of searches are comparable, i.e. ϵ ∼ 10 −3 for invisible limits and ∼ 9 × 10 −4 for visible limits.Therefore, the magenta exclusion band is almost independent of g x .For Fig. 8, where the dark photon has a mass of 10 GeV, the limit from BaBar corresponds to ϵ ∼ 9 × 10 −4 for the dark photon decaying visibly while LEP constrains the invisible decays for ϵ ∼ 4 × 10 −2 .Thus the limits for intermediate g x are rescaled for the right branch solution and become weaker at large values of g x , see Fig. 8 (right) plot.
In Fig. 6 -8, we also show the projected sensitivities of the upcoming accelerator searches with pink dashed lines.For m x = 400 MeV, (which corresponds to m ϕ = 200 MeV here), DUNE is expected to provide the most stringent limits among the near-future visible search experiments [66].In DUNE, a light dark photon is produced in the decays of neutral mesons (π 0 and η).We quote the projected sensitivity of ϵ ∼ 3 × 10 −8 when the dark photon decays 100% into visible particles.On the other hand, LDMX [67], a proposed electron beam dump experiment, is expected to constrain the invisible searches at ϵ ∼ 6 × 10 −5 when x mostly decays invisibly.For a GeV scale dark matter, the best visible search projections are obtained from Belle-II prompt searches into dileptons [66].In the figures, we have used the projected upper limits of ϵ ∼ 2 × 10 −4 for m x = 2 GeV (corresponding to m ϕ = 1 GeV) and ϵ ∼ 1.5 × 10 −4 for m x = 10 GeV (for m ϕ = 5 GeV) from Belle-II.Projections for the invisible searches are taken as ϵ ∼ 3 × 10 −4 for m x = 2 GeV (for m ϕ = 1 GeV) from Belle-II [68] when dark photon decays 100% into invisibles.While realising the future limits, we rescale the numbers according to the branching ratio of dark photon decays, following the same strategy as with the current constraints described before.

Summary plots
The impact of the various constraints on the parameter space of the dark photon model as displayed in Figs. 6 -8 can be summarized as follows.
a. Left branch solutions -Σ 2 0 < Σ2 min For the left branch solutions, large g x values are constrained both by CMB and X-rays while accelerator searches for dark photons constrain the region at large kinetic mixing.When the dark photon is about 400 MeV, a constraint from decays of light mesons into dark photons applies as well, the allowed parameter space corresponds to a narrow region with ϵ roughly between 10 −6 < ϵ < 5 × 10 −4 and g x < 10 −4 .This allowed region shifts toward higher values of g x as the mass of DM increases.In the cases where m ϕ is equal to 1 and 5 GeV, there is an additional region that survives both X-ray and CMB constraints in the ranges 8 × 10 −3 < g x < 3 × 10 −2 (1 GeV) and 2 × 10 −2 < g x < 0.3 (5 GeV).Finally note that the direct detection constraint does not further restrict the parameter space.
b. Right branch solutions -Σ 2 0 > Σ2 min For the right branch solutions which correspond to a larger Σ 2 0 , that is a larger mass gap ∆ = m x − 2 m ϕ , it is the region at low values of ϵ that is constrained by both CMB and X-rays while, as in the previous case, larger values of ϵ are constrained by accelerator searches for dark photons.Direct detection constraints are relevant in the region of parameter space where ϵ and g x are large, especially when m ϕ = 5 GeV.
In the future, direct detection experiments dedicated to detecting low-energy threshold nuclear recoils provide promising limits to constrain DM interaction with SM particles more efficiently.Experiments such as DARK-SPHERE [53]  improvement for a DM mass of 1 (5) GeV, allowing direct detection to significantly probe the allowed region.Again, for a DM mass m ϕ of 1 and 5 GeV, a horizontal band opens up in parameter space, which escapes CMB and X-ray constraints, in the range 8 × 10 −12 < ϵ < 8 × 10 −11 (1 GeV) and 3 × 10 −12 < ϵ < 10 −9 (5 GeV).It is expected that future CMB missions such as PIXIE or PRISM will partially probe these regions.
In the future, various searches for dark photons through their visible or invisible decay modes will be able to probe the region with prompt decays for large values of ϵ and displaced decays for small values of ϵ.Near-future prompt and displaced visible decay search experiments [66] like LHCb, Belle-II, DarkQuest and DUNE are proposed to exclude the kinetic mixing parameter down to ϵ ∼ 3 × 10 −8 for a 200-500 MeV dark photon.Moreover, the sensitivity reach for the invisible decay searches from BaBar is expected to improve by 2 orders of magnitude with LDMX [37] in this mass range.For dark photons mass around the GeV, Belle-II could probe a few 10 −4 with dilepton (visible) [56] and invisible searches [37].

V. CONCLUSION
In this work, we have explored a model where the DM candidate is a GeV-scale scalar species ϕ charged under a new local gauge group U ′ (1).The gauge boson associated to this gauge group, dubbed throughout the article 'dark photon', acts as a vector boson portal between the dark and the SM sectors of the theory.DM annihilation into light SM fermions proceeds through the exchange in the s-channel of this dark photon.In these conditions, DM annihilation has a p-wave behavior at low energy.However, it can be significantly enhanced if two conditions are met.First, the dark photon mass m x must be slightly larger than twice the DM mass m ϕ and, second, the decay width of the dark photon must be smaller than twice the mass gap ∆ between m x and 2m ϕ .If both conditions are met, a Breit-Wigner resonance appears at a DM dispersion velocity directly related to this mass gap, i.e. when the DM temperature is equal to ∆.The smaller the mass gap, the smaller the DM dispersion velocity at peak annihilation.
This has profound implications for the cosmological behavior of our DM candidate.In most models in the literature, the bulk of DM annihilation takes place during the freeze-out process.Shortly after freeze-out, DM annihilation stops and the cosmological abundance of DM reaches rapidly its relic level.In our model, on the contrary, the peak of DM annihilation may occur well after freeze-out if the mass gap is very small.Moreover, the relic abundance is always shaped at that particular moment.Determining the evolution of the DM temperature during the expansion of the Universe turns out to be crucial, insofar as DM annihilation mostly occurs when that temperature is of order the mass gap ∆.That is why kinetic decoupling between DM and the primordial plasma must be dealt with.As we showed, failure to do so results in a DM relic abundance being overestimated by several orders of magnitude.
At fixed DM mass m ϕ , dark charge g x and mixing angle ϵ, we found in general two values for the dark photon mass m x for which the DM relic abundance is equal to the cosmological measurement [1].The smallest value, for which Σ 2 0 < Σ2 min , is classified as the left branch solution while the largest value, for which Σ 2 0 > Σ2 min , is classified as the right branch solution.For both possibilities, we explored the (g x , ϵ) plane for a DM mass m ϕ set equal to 200 MeV, 1 GeV and 5 GeV.For small values of g x and ϵ, the scalar DM candidate overshoots the cosmological observed value.For left branch solutions only, large values of g x and ϵ are also excluded since in this case scalar DM undershoots the measurement.We then applied several constraints to the surviving regions.In the domains where ϵ is small, the peak of DM annihilation occurs when the CMB energy spectrum is most sensitive to energy injection.We found µ and y-distortions exceeding by far the bounds set by observations.For smaller values of g x (left branch) and larger values of ϵ (right branch), the regions are excluded by X-ray observations.These astrophysical constraints are complemented by the limits set by colliders and direct detection experiments for large values of ϵ.
In all the cases which we investigated, we found regions having successfully passed all the tests.In the left branch case, this region lies at the upper-left corner of the band where g x is small and ϵ is large.In the right branch case, the surviving domain corresponds roughly to values of ϵ in the range between 10 −7 and 10 −3 .In these allowed regions, the mass gap ∆ is not too small and the peak of DM annihilation occurs shortly after freeze-out.
Although energy injection at late cosmological times is severely constrained, we found that new domains open up in parameter space, all the more so when DM is heavy.Vertical (horizontal) bands appear in the left (right) branch plots of Figs. 7 and 8.We expect the future CMB missions PIXIE (NASA) or PRISM (ESA) to partially close these possibilities since the future instruments will reach a sensitivity of 10 −8 on both y and µ-distortions [69].We anticipate that these regions will nevertheless survive.
In the future, the region at large values of g x and ϵ could be further probed with dedicated low-energy threshold direct detection experiments such as DARKSPHERE [53] and SBC [54,55] and improve the current limits by a few orders in our mass range of interest.Moreover, several planned searches for dark photons through their visible and invisible decays at accelerators such as LHCb, Belle-II, DUNE, DarkQuest or LDMX [56,66] are proposed to improve the upper bound on ϵ by a few orders of magnitude for a few hundreds of MeV and by an order of magnitude for a few GeV dark photon.We have discussed this in details in the text.
Our results should finally be taken with a grain of salt.Our conclusions are based on the key assumption that DM is always in thermal equilibrium with itself.In the regions having passed all our tests, kinetic decoupling occurs well after DM has become non-relativistic.Actually, before kinetic decoupling happens, DM reaches a state of inner thermal equilibrium through its close contact with the primordial plasma, exactly like photons do with electrons just before recombination.But, after kinetic decoupling, our assumption that DM is still in thermal contact with itself needs to be scrutinized.Such an investigation was beyond the scope of this exploratory analysis.In a follow-up study, we plan to investigate this point by solving, for instance, the transport Boltzmann equation for DM species and evolve in time the particle momentum distribution [70].A more ambitious goal would be to embed the Lagrangian II.1 in a more general set-up and to study the cosmological consequences of the overall theory.This would allow to determine whether or not the DM scalar ϕ is actually in thermodynamical equilibrium with the primordial plasma before it becomes non-relativistic.Beyond that point, J drops exponentially until the third regime is reached.(iii) At large values of |a|, the integral J decreases actually as 1/a 2 .In the right panel, the curves for which the ratio is small, exhibit a Breit-Wigner enhancement, while for larger values of the ratio, we observe the smooth transition from the constant value of 1 to the asymptotic form.
The averaged cross-section ⟨σ ann v⟩ scales this time as Σ 2 .The annihilation is p-wave suppressed at low velocities.
• a < 0 : In this case, the dark photon is heavier than a ϕ-φ pair at rest.The annihilation can be enormously enhanced if the velocities of the incoming DM scalars can fill the gap ∆ between 2m ϕ and m x and if the decay width of the dark photon that mediates the process is very narrow.The annihilation proceeds then through a Breit-Wigner resonance.The smaller is Λ 0 with respect to Σ 0 , the narrower and higher the resonance.The peak is actually reached when the one-dimensional dispersion velocity Σ of the DM scalars is of order Σ 0 .
We remark with reference to Eq. II.11 that when Σ increases, both |a| and b decrease while the ratio |a|/b remains fixed at Σ 2 0 /Λ 2 0 .The behavior of the integral J is now more involved than when a was positive.The left panel of Fig. 9 features the evolution of J as a function of |a| for a fixed ratio b/|a| of 10 −4 .Three regimes are clearly visible : 1.When the parameter |a| goes to 0, the integral converges toward 1 as it should, since J(0, 0) = 1. 2. For small values of the ratio b/|a|, and as long as |a| is not too large with respect to 1, the Breit-Wigner resonance sets in.The region of integration which contributes most to J is centered around t = |a|.It corresponds to the annihilation peak whose full width at half-maximum is 2b ≪ |a|.We can write the integral as where the variable u = t − |a|.The integral J is well approximated by the function The evolution of J as a function of |a| is plotted in the right panel of Fig. 9 for several values of the ratio b/|a|.When this ratio is small, i.e. when the Breit-Wigner resonance through which the annihilation proceeds is narrow, this translates into Λ 0 being much smaller than Σ 0 and the integral J has the same behavior as in Fig. 9.It features a bump whose highest point is reached at |a| = 1/2.On the contrary, when Λ 0 is larger than Σ 0 , i.e. for large values of the ratio b/|a|, the contribution from J 1 becomes subdominant and the integral J does not exhibit any enhancement.It evolves smoothly from the constant value of 1 to its asymptotic form J 2 .The transition between these regimes takes place for |a| ∼ |a|/b.We recover the same behavior as if a were positive.

Appendix B: Scalar Dark Matter thermalization
In this section, we present a simplified analysis of the thermalization of scalar DM with the SM plasma.A complete treatment would require the knowledge of the overall theory, which is beyond the scope of our exploratory work.We will concentrate here on the interactions between scalar DM and SM fermions that are encoded in Lagrangian II.1, the starting point of our analysis.
We will also assume that scalar DM is thermalized with itself, and that a DM temperature T ϕ can be defined at all times.This is certainly true whenever collisions between DM scalars and SM fermions are rapid enough to ensure efficient energy exchange between these populations and to establish thermalization.We will assume that DM reaches a state of inner thermal equilibrium after kinetic decoupling from the SM plasma has occured, allowing the DM temperature T ϕ to be defined also in this situation.This may appear as an oversimplification.Going beyond it would make the problem orders of magnitude more complicated.We would have to solve directly the Boltzmann equation and study the evolution in time of the DM distribution function in momentum space [70].We will defer such an investigation to a future work.Assuming that DM reaches inner thermal equilibrium with temperature T ϕ leads already to a particularly rich and complex phenomenology, which could be the starting point for further investigations.
The question of DM thermalizing with the primordial plasma is paramount insofar as the annihilation crosssection II.10 crucially depends on the DM dispersion velocity Σ and DM temperature T ϕ ≡ m ϕ Σ 2 .As encrypted in Lagrangian II.1, scalar DM exchanges energy with the SM plasma through annihilations into, recreations from and collisions upon SM light fermions.Our aim is to model these processes to establish an equation that drives the evolution of the DM temperature T ϕ with respect to the SM plasma temperature T .Notice that the DM heat capacity is small compared to that of the SM plasma.As it becomes non-relativistic, DM annihilates and its density actually drops, hence a negligible contribution to the overall heat capacity of the primordial plasma.The temperature of the latter still evolves at constant entropy, decreasing approximately like a −1 , where a is the scale factor of the Universe.The SM plasma is not affected by its thermal contact with DM, or the breaking of it.This is not so for the DM temperature which relaxes rapidly toward the SM temperature at early times.At some point, called thermal or kinetic decoupling, this relaxation slows down, and T ϕ cannot follow T anymore.The DM temperature decreases afterward more rapidly than the plasma temperature.
In Sec.B 1, we show that DM can be thermalized through its annihilation into and recreation from SM fermions.We show that kinetic decoupling occurs always after freeze-out.In Sec.B 2, we model the energy exchanged between DM and the SM plasma through collisions.We eventually establish the master equation for T ϕ in Sec.B 3 and explain how we find the kinetic decoupling point.

Thermalization through annihilations
We investigate here whether or not the thermalization of ϕ and φ particles with the SM plasma could be ensured solely through their annihilation into, and production from, standard model fermions We want to determine if DM species are thermalized, and thermal contact is established with the SM plasma, should this reaction be fast enough.We assume that DM particles are thermalized with each other so that a DM temperature T ϕ can be defined.We would like to determine how fast T ϕ relaxes toward the plasma temperature T .For this, let us consider a volume V of space.It is filled with DM species ϕ and φ whose densities n ϕ and n φ are equal insofar as no asymmetry is assumed.Let us define n ≡ n ϕ ≡ n φ.The number of DM particles inside the volume V is The DM energy stored inside that volume is We assume DM to be non-relativistic as we are dealing with the freeze-out process at temperatures below m ϕ .The pressure of the DM gas is P = (n ϕ + n φ)T ϕ = 2nT ϕ .During the lapse of time dt, the number N of DM species inside volume V varies according to the chemical reaction B.1 Two DM species disappear or are created per reaction.The creation term can be derived by noticing that it should cancel the annihilation term at thermodynamical equilibrium.The density n eq corresponds to a population at temperature T with vanishing chemical potential.In the non-relativistic regime under consideration, it is given by Eq.III.9.
On the other hand, relation B. At early times, before freeze-out starts, we can asume DM to be in thermodynamical equilibrium with the rest of the plasma so that n = n eq while T ϕ = T .The annihilation rate Γ ann is given by ⟨σ ann v⟩ T n eq and is equal to the rate Γ F As both Eqs.B.13 and B.14 are similar, we conclude that freeze-out (aka chemical decoupling) is concomitant with thermal decoupling (aka kinetic decoupling).We can even guess that freeze-out occurs slightly before kinetic decoupling.Actually, the codensity N relaxes toward its dynamical equilibrium value N eq with rate Γ ann , while the DM temperature T ϕ relaxes toward the plasma temperature T with the slightly larger rate Γ ann + 2H.We also notice that the right hand side term Γ ann N eq of the codensity equation B.14 drops much faster than Γ ann T , its temperature counterpart in equation B. 13.The codensity N eq has an exponential dependence exp(−m ϕ /T ) and decreases much faster than T .
To conclude, we have proved that DM annihilation into, and production from, SM light fermions is able alone to thermalize DM with the plasma.Kinetic decoupling in that case occurs slightly after freeze-out.Reaction B.1 results from the interaction Lagrangian II.1 which also implies the existence of collisions between DM and light fermions discussed in the next section.

Thermalization through collisions with SM species
There are of course collisions between light SM fermions and DM.These contribute to the thermalization of DM since energy is exchanged between both populations.Let us develop a simplified calculation of the energy transferred from the fermions f to the DM species.The latter are non-relativistic since we are interested here in the period of kinetic decoupling which occurs after freeze-out, for a plasma temperature below m ϕ .We can safely treat the DM scalars as if they were at rest in the cosmological frame and compute their recoil energy as they are impacted by incident fermions.We focus on the collision f + ϕ −→ f + ϕ . (B.15) A dark photon is exchanged in the t-channel between the fermionic line and the scalar line.The fermion f is ultrarelativistic as we are interested this time in plasma temperatures T above the mass m f .In the opposite situation, pair annihilation drives the population of fermions f to extinction, with number density n f dropping like exp(−m f /T ), and energy transfer stops.The incident fermion in reaction B.15 has energy ϵ f and is scattered through the angle θ with respect to its initial direction.The scalar ϕ, initially at rest, recoils with kinetic energy E R .A straightforward but exact calculation yields We recover the well-known relation of the Compton effect.In the regime where the plasma temperature T is significantly smaller than the DM mass m ϕ , so is the average fermion energy.The parameter u is small with respect to 1 and the recoil energy boils down to The collisions of the fermions f and f on the scalars ϕ and φ result in the increase of the DM internal energy contained in volume V with rate where E R (Ω) is given by B.17.As fermions are ultra-relativistic, their velocity v f is equal to the speed of light.Their number density obeys Fermi-Dirac statistics  We can apply detailed balance to describe now the energy flowing from scalar DM to fermions as a result of collisions.This boils down to replacing T 7 by T 7 ϕ in the previous expression.Going a step further by linearizing the net gain of DM internal energy per unit time with respect to the temperature difference T ϕ − T , we get The coefficient C col is given by where We have summed over all fermionic populations which are ultra-relativistic when the plasma temperature is T , hence the effective charge Q 2 eff .Notice that we have so far concentrated on fermions but we might be in the situation where charged pions also collide on DM scalars.This may happen at the end of the quark-hadron phase transition, although pions are not strictly ultra-relativistic at that time.We can slightly modify the definition of Q 2 eff to include these charged scalar species.Their collision cross-section is enhanced by a factor of 3 with respect to fermions and the statistical factor must also be modified.We describe now Q 2 eff as the smooth function of the plasma temperature T where g π = 2 and Q 2 π = 1.

Kinetic decoupling
Taking into account both annihilations and collisions, the evolution of the scalar DM temperature becomes At early times, Γ KD rel exceeds Γ KD eq by orders of magnitude and DM is thermally connected to the SM plasma.But as time goes on, Γ ann drops exponentially like exp(−m ϕ /T ) while Γ col drops like T 6 .The relaxation rate Γ KD rel decreases faster than the rate Γ KD eq at which the kinetic equilibrium evolves.Kinetic decoupling occurs when both rates are equal Γ KD rel (x KD ) = Γ KD eq (x KD ) , (B.31 where the kinetic-decoupling point x KD is defined as the ratio m ϕ /T KD .Afterward, the DM temperature follows the equation and T ϕ decreases like a −2 , where a is the scale factor of the Universe.DM behaves then like a non-relativistic gas undergoing adiabatic expansion.

FIG. 1 :
FIG.1:The variation of the annihilation cross-section as a function of the dispersion velocity Σ for different values of b/|a| = Λ 2 0 /Σ 2 0 .On the horizontal axis, the rescaled variable (2m ϕ /mx)(Σ/Σ0) has been used.When Λ0 is smaller than Σ0, the cross-section is enhanced by a Breit-Wigner resonance.Above a velocity of order Σ0, where its peak value is reached, ⟨σannv⟩ drops like Σ −3 to reach the asymptotic behavior Σ −2 .Below the peak, the p-wave annihilation regime sets in and ⟨σannv⟩ is proportional to Σ 2 .For large values of Λ0 with respect to Σ0, the two asymptotic regimes only appear.

FIG. 4 :
FIG.4: Allowed parameter space in the (gx, ϵ) plane, for left branch solutions (Σ 2 0 < Σ2 min ), and for a DM mass of 1 GeV.From the upper-left to the lower-right panel, the color code shows xKD, Σ 2 0 , the corresponding DM annihilation cross-section in the Milky Way ⟨σannv⟩MW and the CMB µ-distortion.In the bottom-left panel, the hatched region features the exclusion constraints by XMM X-ray measurements.In the bottom-right panel, the forbidden hatched region is drawn from CMB µ-distortion as observed by FIRAS (COBE).It encompasses the domains excluded by CMB y-distortions (dashed line) and anisotropies (dotted line).The full interpretation of these plots is provided in the text.

FIG. 6 :
FIG. 6: Summary plot for m ϕ = 200 MeV: constraints on the parameter space (gx, ϵ).The left and right panels correspond to the left (Σ 2 0 < Σ2 min ) and right (Σ 2 0 > Σ2 min ) branch solutions, respectively.The color scale indicates the resulting ⟨σannv⟩MW.The colored patches show excluded regions from CMB µ-distortion constraints (blue), X-rays measured by XMM (gray), colliders (pink) and direct detection (orange).Projected exclusion limits are depicted for future direct detection experiments (orange dashed lines) and future accelerator searches (pink dashed lines).More details are provided in the Sect.IV.

FIG. 9 :
FIG. 9: The integral J(a, b) is plotted as a function of |a| (a < 0) for a fixed ratio b/|a| of 10 −4 (left) and several values of the ratio b/|a| (right).In the left panel, three behaviors are clearly visible -(i) When the parameter |a| goes to zero, the integral J tends to 1. (ii) The Breit-Wigner regime sits in the intermediate region where J steadily increases like |a| until |a| reaches 1/2.Beyond that point, J drops exponentially until the third regime is reached.(iii) At large values of |a|, the integral J decreases actually as 1/a 2 .In the right panel, the curves for which the ratio is small, exhibit a Breit-Wigner enhancement, while for larger values of the ratio, we observe the smooth transition from the constant value of 1 to the asymptotic form.

9 ) 10 )
where H = ȧ/a is the expansion rate and a is the scale factor of the Universe at time t.Combining both equalities B.4 and B.5 yields the well-known equation dn dt = −3Hn − ⟨σ ann v⟩ T ϕ n 2 + ⟨σ ann v⟩ T n 2 eq .(B.6)Alternatively, we can recast relation B.4 into dN dt + ⟨σ ann v⟩ T ϕ n N = {⟨σ ann v⟩ T n eq } N eq .(B.7)As long as the annihilation rate Γ ann ≡ ⟨σ ann v⟩ T ϕ n (B.8) is larger than the rate with which the right hand side term evolves, a dynamical equilibrium is reached and Ṅ can be disregared in equation B.7.If both temperatures T ϕ and T are furthermore equal, the DM density n is given by its chemical equilibrium value III.9.During the lapse of time dt, the DM internal energy B.3 varies by an amount dU = −P dV − 2 ⟨σ ann v⟩ T ϕ n 2 V dt m ϕ + 3 2 T ϕ + 2 ⟨σ ann v⟩ T n 2 eq V dt m ϕ We have also applied detailed balance, noticing that each time a pair of DM scalars disappears, an average energy twice equal to m ϕ + 3/2 T ϕ is removed from the DM population.The amount of energy given to DM each time a pair of scalars ϕ φ is created is twice equal to m ϕ + 3/2 T , with this time T instead of T ϕ .Actually when thermodynamical equilibrium is reached between DM and the SM plasma, the only change in the DM internal energy comes from the pressure work −P dV.The variation of the DM internal energy B.3 can also be written as dU = dN m ϕ Using expressions B.7, B.9 and B.10, we derive the equation fulfilled by the DM temperature dT ϕ dt = −2HT ϕ − ⟨σ ann v⟩ T n 2 eq n (T ϕ − T ) .(B.11) rel defined in Sec III A. The DM temperature evolves asdT ϕ dt = −2HT ϕ − Γ ann (T ϕ − T ) .(B.12)This equation can be recast into dT ϕ dt + {Γ ann + 2H} T ϕ = Γ ann T , (B.13) with the same structure as B.7 which boils down to dN dt + Γ ann N = Γ ann N eq .(B.14)

2 x ϵ 2 e 2 m 4 x m ϕ T 6 .
dT ϕ dt = −2HT ϕ − ⟨σ ann v⟩ T n 2 eq n + C col T 6 (T ϕ − T ) .(B.26)For numerical resolution purposes, this equation can be recast into a form similar to relation III.4dT ϕ dt + {Γ ann + Γ col + 2H} T ϕ = {Γ ann + Γ col } T , (B.27)whereΓ ann = ⟨σ ann v⟩ T n 2 eq n ≃ ⟨σ ann v⟩ T n eq = Γ F rel while Γ col = C col T 6 = A col Q 2 eff g (B.28)The DM temperature T ϕ relaxes toward the plasma temperature T with rateΓ KD rel = Γ ann + Γ col + 2H , (B.29)while the rate of evolution of kinetic equilibrium can be defined as in III.6 through the identity Γ KD eq ≡ d dt ln {(Γ ann + Γ col ) T } .(B.30) 1 , (B.19) where g f denotes the spin degeneracy factor.The Feynman diagram associated to reaction B.15 yields the differential collision cross-section Plugging the two last relations into the rate of collisional energy transfer B.18 yields