Indirect searches for dark matter bound state formation and level transitions

Indirect searches for dark matter (DM) have conventionally been applied to the products of DM annihilation or decay. If DM couples to light force carriers, however, it can be captured into bound states via dissipation of energy that may yield detectable signals. We extend the indirect searches to DM bound state formation and transitions between bound levels, and constrain the emission of unstable dark photons. Our results significantly refine the predicted signal flux that could be observed in experiments. As a concrete example, we use Fermi-LAT dwarf spheroidal observations to obtain constraints in terms of the dark photon mass and energy which we use to search for the formation of stable or unstable bound states.


Introduction
Most of the dark matter (DM) research in the past decades has focused on DM with contact-type interactions, i.e. interactions mediated by particles of similar or larger mass than the DM itself, m med m DM . Indeed, in the prototypical WIMP (Weakly Interacting Massive Particle) scenario, DM was envisioned to couple to the weak interactions of the Standard Model (SM) and have mass m DM ∼ m W,Z ∼ 100 GeV. The current collider, direct detection, and indirect detection searches strongly constrain this scenario. Nevertheless, they still allow for WIMP DM around or beyond the TeV scale. The same conclusion essentially holds for a variety of models in which DM communicates with the SM particles via non-SM mediators. However, for WIMP DM with m DM TeV m W,Z , the weak interactions manifest as long-range [1].
On a more fundamental and model-independent level, the unitarity of the S-matrix suggests that the long-range character of the interactions is a generic feature of viable thermal-relic DM models in the multi-TeV mass range and above [2]. Indeed, unitarity sets an upper bound on the partial-wave inelastic cross sections, whose physical significance is the saturation of the probability for inelastic scattering. This in turn implies an upper bound on the mass of DM produced via thermal freezeout [3], of the order of 100 TeV [2,4]. The parametric dependence of the unitarity limit on the inelastic cross sections shows that the 100 TeV regime can be approached or reached only by interactions that manifest as long-range [2].
Long-range interactions imply the emergence of non-perturbative effects that can affect significantly the DM phenomenology. A long-range force distorts the wave function of a pair of DM particles, and consequently affects all their interaction rates at low velocities. This is the well known Sommerfeld effect [5,6], which has been extensively studied in the DM literature, both in WIMP and hidden-sector models (see e.g. [1,2,[7][8][9][10][11][12][13][14][15][16][17][18][19][20]). It has been shown to decrease the DM density in the early Universe, and enhance the expected indirect detection signals. Another potentially more consequential implication of long-range interactions is the existence of bound states [4,9,21]. Bound-state formation (BSF) -which is also affected by the Sommerfeld effect -can alter the DM phenomenology in a variety of ways.
To delineate the consequences of DM bound states, we discern two broad categories.
Unstable bound states. The formation of particle-antiparticle (positronium-like) bound states that can decay into radiation opens a new two-step DM annihilation channel. In models that feature co-annihilation between different species and/or non-Abelian forces, a variety of unstable bound states may exist. The formation and decay of unstable bound states diminish the DM density in the early Universe [4], thereby altering the expected DM mass and couplings and affecting all experimental signatures [4,[22][23][24][25][26][27][28][29][30]. During the CMB (Cosmic Microwave Background) period and inside galaxies today, the bound state decay products enhance the high-energy radiative signals searched by telescopes [9,18,20,[31][32][33][34][35].
Stable bound states. The formation of stable bound states alters -typically screens or curtailsthe DM self-interactions inside halos [36,37], which are expected to affect the galactic structure [38,39]. Moreover, stable bound states affect the DM direct detection signatures [40,41]. Stable bound states arise typically either due to confining forces (hadronic bound states), and/or due to weak forces in models of asymmetric DM. The latter scenario hypothesizes that the DM relic density is, analogously to ordinary matter, due to an excess of dark particles over antiparticles [42][43][44] that cannot be annihilated in the early Universe even if the DM annihilation cross section is very large. It follows that DM may possess significant couplings to light force mediators that in turn may render BSF rather efficient.
The DM capture into bound states, be they stable or unstable, invariably necessitates the dissipation of energy. The amount dissipated -of the order of the binding energy of the bound stateis significantly lower than that radiated in the typical DM annihilation and decay of unstable bound states. Pearce and Kusenko first suggested that the energy dissipated during BSF at late times may give rise to novel signals detectable via indirect searches [45]. Transitions between bound levels may also produce similar signals. Indeed, the available energy may be dissipated radiatively, most commonly via emission of a force mediator. In multi-TeV WIMP models, BSF inside halos can occur via emission of a mono-energetic photon in the multi-GeV range [35]. In models where DM couples to non-SM forces, the emitted mediators may decay into SM particles whose cascades produce a more extended spectrum of photons and other stable SM particles [18,20,46]. 1 The purpose of the present work is to initiate a systematic investigation of indirect constraints on BSF.
The indirect signals emanating from BSF and other level transitions can provide a powerful probe of asymmetric DM models, where late-time DM annihilation is highly suppressed due to the absence of antiparticles [20,[50][51][52], unless a mechanism exists that erases the asymmetry at late times [53][54][55]. In contrast to annihilation, and because asymmetric DM can accommodate large DM-mediator couplings, BSF can be quite efficient. Part of the parameter space where this occurs is in fact interesting for an additional reason, that it provides a viable framework of self-interacting DM [20,36]. Indirect signals from the formation of stable bound states in asymmetric DM models have been proposed in Refs. [45,46,[56][57][58]. As a concrete example that we consider below in more detail, we mention here the formation of dark atoms via emission of light dark photons that subsequently decay into SM particles via kinetic mixing with hypercharge [46]. While the parameter space of the model is broader, it has been shown that dark atomic transitions between levels with MeV-scale splittings could inject low-energy positrons in the Milky Way at a sufficient rate to account for the observed 511 keV line [46].
Even in the context of symmetric or self-conjugate DM, BSF signals may provide an important probe, since they may exhibit different spectral features and resonant structure than direct annihilation [59]. Moreover, for very heavy DM whose annihilation signals fall outside the energy range of the various telescopes, the low-energy radiation could fall within the energy range probed by telescopes, and could thus be employed to constrain a wider range of DM masses. 2 The radiative BSF cross sections can be comparable to or even significantly larger than the direct annihilation cross sections [4,21,23,24,29,59]. In fact, the BSF cross sections in galactic environments may exceed the so-called canonical annihilation cross section, σv rel ≈ 3 × 10 −26 cm 3 /s, by orders of magnitude due to different reasons. These include a large Sommerfeld enhancement at low velocities, and possibly the associated parametric resonances in the case of massive mediators [59], as well as, in the case of asymmetric DM, a larger DM-mediator coupling than that required to attain the observed DM density via freeze-out in the symmetric limit [2]. However, the accurate estimation of the expected BSF signals, and indeed of any DM experimental signature, necessitates computing the cosmology first [43]. If bound states exist, then they may form efficiently in the early Universe. As already mentioned, the formation and decay of unstable bound states in the early Universe decreases the DM density [4], and therefore alters the predicted DM parameters that determine the late-time BSF rate [18]. The formation of stable bound states in the early Universe changes the density of particles available to participate in the corresponding processes inside halos, and thus again affects the expected indirect signals [46].
The structure of the paper is as follows. To flesh out the above, we begin in Section 2 by introducing an atomic DM scenario with a light albeit massive dark photon that mixes kinetically with hypercharge. After summarising the cosmology of the model, we estimate the indirect signals expected from the DM capture into dark atoms via emission of dark photons. Compared to previous studies [46], we employ improved numerical calculations of the BSF cross sections [59], and compute the γ-ray flux from the cascades of the charged particles produced in the dark photon decays. In Section 3, we briefly consider the recasting of existing constraints on DM annihilation into SM particles for the purpose of constraining BSF, before deriving new constraints on BSF occurring via dark photon emission. We use Fermi-LAT observations of dwarf spheroidal (dSph) galaxies, which provide a DM rich environment with relatively lower background compared to the Galactic Centre. The constraints are cast in terms of the DM mass and the energy dissipated, such that they can be used in models with different underlying dynamics. They are applicable to BSF, as well as excitation processes occuring via DM collisions and followed by de-excitations. The predictions of the atomic DM model of Section 2 are confronted with the derived constraints in Section 4, where we also discuss p e V H p e H V · · · · · · Figure 1: The process targeted in this project. A dark proton and dark electron combine to form dark hydrogen through the emission of a massive mediator (dark photon). The mediator decays into SM particles through its kinetic mixing with hypercharge, which eventually yields lower energy photons, which can be searched for by Fermi.
further applications. Some general remarks are then drawn in the conclusion.
2 Atomic dark matter with a massive dark photon

The model
We assume that DM is charged under a dark U (1) D gauge symmetry, and that it carries a particleantiparticle asymmetry conserved at low energies due to a global dark baryonic symmetry governing the interactions of the dark sector. If U (1) D is unbroken, then gauge invariance mandates that there must be at least two dark particle species with compensating asymmetries, such that the dark electric charge of the Universe vanishes. This remains true if the dark photon acquires a mass via the Stückelberg mechanism, or if U (1) D is broken via a Higgs mechanism that operated in the early Universe after the dark baryogenesis took place. The latter implies that the generated dark photon mass is sufficiently small. We refer to [43] for the detailed considerations.
Considering the above, we will assume that the dark photon V has a small non-zero mass m V , and that DM consists of two species of fundamental Dirac fermions, the dark protons p and the dark electrons e, with opposite charges and masses m p m e . Thus, the low-energy physics of the dark sector we explore in this paper is summarised by the Lagrangian: The covariant derivative is D µ = ∂ µ ± ig D V µ for p and e respectively. The field strength tensor is is the dark fine-structure constant. The dark photons may decay into SM particles through the kinetic mixing with hypercharge, controlled by the dimensionless parameter . Here c w ≡ cos θ w where θ w is the Weinberg angle. Constraints on the dark photon and on DM direct detection via dark photon exchange are compiled in Appendices A and B respectively.
High energy completions of this scenario, including mechanisms for the generation of the dark matter-antimatter asymmetry, that could be potentially related to that of ordinary matter, can be found e.g. in Refs. [63][64][65][66][67], and the DM freeze-out has been previously studied in Refs. [2,20]. Here we shall only use that the dark proton-antiproton and dark electron-positron asymmetries are equal, and that the dark antiparticles were efficiently annihilated in the early Universe with an equal amount of dark particles, thereby leaving a Universe that contains globally (nearly) equal densities of dark protons and dark electrons, n p ∼ = n e ≫ np, nē. The exact number of residual antiparticles depends on the effective annihilation cross section in the early Universe, here controlled by the coupling α D ; in order for the DM density to be set largely by the primordial asymmetry, it is sufficient that α D is somewhat higher than that for symmetric thermal relic DM of the same mass [2,20,50].
The symmetric DM realisation of this scenario, containing only one dark species, has been studied in Refs. [18,68], with particular emphasis on the indirect constraints due to late-time DM annihilations. Remarkably, annihilation constraints arising from the small but non-zero residual density of dark antiparticles, exist also in the asymmetric regime for late-time asymmetries np/n p 10 −3 , due to the large Sommerfeld enhancement of the annihilation cross section that compensates in part for the suppression of the annihilation rate due to the small residual density of antiparticles [2,20]. For larger asymmetries, i.e. larger values of α D , the annihilation rate falls below the sensitivity of the current observations. Nevertheless, a larger α D implies that the formation of stable bound states may be possible for a larger range of m p , m e , giving rise to radiative signals [46] that we shall now explore.
If the dark photons are sufficiently light, then the dark protons and the dark electrons can form dark hydrogen atoms. The capture into atomic bound states may occur via emission of a dark photon, provided it is kinematically allowed, and as illustrated in Fig. 1. The dark atom formation may occur in the early Universe (dark recombination), as well as at late times, during the CMB period or inside galaxies today. In the following, we specify the relevant BSF cross section, the ionized fraction of DM that may participate in this process today, the branching fractions of the dark photons into SM particles, and the γ-ray spectrum resulting from the cascades of the latter.

Formation of dark atoms
For convenience, we define the following parameters (2.5) The first is important in determining the strength of the Sommerfeld enhancement and the overlap of the scattering-state and bound-state wave functions. The second is the ratio of the dark atom Bohr radius to the range of the dark-photon-mediated interaction, and parametrises how long range this interaction manifests. The last parameter is the p − e reduced mass. Bound levels of pe pairs exist if ξ > 0.84 [59]. They may form radiatively, via emission of a dark photon. We will consider capture into the ground-state only, which is the dominant BSF process and most exothermic transition, with the energy available to be dissipated being [21] where E D ≡ γ 2 D (ξ) × µ D α 2 D /2 is the absolute value of the binding energy. The factor γ D (ξ) 1 parametrises the departure from the Coulomb value. The cross section for capture into the ground state is [59] (σv rel ) where parametrises the phase-space suppression due to the massive dark photon. Both γ D (ξ) and the function S BSF (ζ, ξ) in Eq. (2.7) are computed numerically according to Ref. [59]. In the Coulomb limit ξ → ∞, the latter takes the analytical form [21,36,59] S BSF (ζ) 2πζ 1 − e −2πζ 2 10 3 (2.9) In fact, the Coulomb approximation is satisfactory for µ D v rel m V , or equivalently ξ ζ [59]. We note that the capture into the ground state is a p-wave process. While in the Coulomb regime and for α D v rel , (σv rel ) BSF exhibits the characteristic Sommerfeld scaling ∝ 1/v rel as seen from Eq. (2.9), the finite m V implies that at velocities v rel m V /µ D , the BSF cross section recovers the perturbative scaling (σv rel ) BSF ∝ v 2 rel [59]. An example of the cross section including the effects of the finite mediator mass is shown in Fig. 2. Because the BSF cross section is suppressed at v rel α D , as seen from Eq. (2.9), the phase-space suppression (2.8) implies that ξ 1 whenever BSF is kinematically allowed and significant, thus to a good approximation γ D (ξ) 1. At threshold, the dark photon is produced at rest and the ratio of b T : b L reads 2 3 : 1 3 . In the limit of E D + µ D v 2 rel /2 m V we instead recover 1 : 0, as anticipated since the longitudinal polarisation is unphysical in the limit m V → 0. The significance of this is that the angular distribution of dark photon decay products and the final-state ordinary photons depend on the dark-photon polarisation. Hence, once we boost from the dark photon rest frame into the observer frame, the energy spectrum of the final-state photons will also depend on the polarisation. The branching ratios (2.10) are illustrated in Fig. 3 for some favourable choices of parameters for returning a significant indirect detection signal.

Residual ionisation
The formation of dark atoms may also occur in the early Universe, thereby reducing the density of dark ions today. To compute the BSF rate inside galaxies, we must therefore consider the ionized fraction of DM. This is defined as f ion ≡ n p n H + n p , (2.11) where n p is the number density of unbound dark protons and n H is the number density of dark hydrogen. After dark recombination in the early Universe, the residual ionized fraction can be estimated under the assumption of Saha equilibrium and freeze-out as where we have included the phase space factor of the BSF cross section (2.7). The factor τ rec = min[1, T D /T SM ] rec takes into account the potentially different temperatures of the dark sector and the SM plasma during dark recombination (cf. Appendix C), which occurs at T D ∼ 0.007E D [69]. A more detailed computation of the dark recombination that takes into account multi-level transitions has been performed in Ref. [69], according to which the approximation of Eq.  that dark atom-atom or atom-ion collisions inside halos partially reionise DM. If this is significant (cf. Footnote 4), then the use of Eq. (2.12) underestimates the indirect signals due to BSF, thus leading to conservative constraints. Considering the above, and given also the various other sources of error -notably the J-factors -that will enter our analysis below, we shall proceed using Eq. (2.12) throughout. An example of f ion is shown in Fig. 4.

Interactions between the species inside halos
If p − e interactions are sufficiently strong, then the dark electrons may receive a kick and escape the DM halo. This would in turn suppress the expected BSF rate and indirect signals. One may wonder if the build-up of a net charge in an area of the halo is consistent with the long-range Yukawa interaction. However, for the mediator masses considered here, m V > MeV, the range of the interaction is at most of the order of picometer. Hence ejected electrons will not be drawn back into the region by the dark gauge force.
The interaction between the DM species can be found from the formulas given in [70], which take into account the possible long-range interaction due to the light mediator, by making the replacement m X /2 → µ D . Using these cross sections, we estimate the p − e scattering rate in the areas of the halo of interest. For fully ionised DM, the scattering rate of a dark electron on dark protons is given by We take the typical DM density of a dSph at the region of interest to be ρ DM ≈ 10 GeV/cm 3 [71], 3 and the typical relative velocity v rel ∼ 20 km/s [72]. Assuming a lifetime of 10 billion years, in Fig. 5 we show parameter regions where the electrons undergo on average one or more scatterings and could therefore thermalise. We see that for m V O(0.01) GeV, the thermalisation is inefficient and we expect that the density of dark electrons in the halo is essentially the same as that of the dark protons. This rough bound on m V is not affected much by considering a larger ρ DM in Eq. (2.13) because the elastic cross section drops very rapidly with increasing m V . Note the estimate can be refined by including f ion in Eq. (2.13), which would reduce Γ scat , although we should then also take into account the H − e scatterings, whose cross section is however more suppressed due to screening. Considering our later results in Section 3, and the pre-existing constraints on dark photons summarised in Appendix A -which allow mostly for m V O(0.1) GeV -it becomes clear that inclusion of such effects will not change the parameter space of interest in the present study. We therefore do not include such complications here. 4

Dark photon decay
To obtain the expected signal we need the dark photon branching fractions. A standard perturbative, tree-level, calculation allows one to find the partial widths into the individual decay channels using the couplings of the dark photon to the SM fermions (2.14) The couplings to the chiral components of the fields are given by [73] c where Here c α ≡ cos α, s α ≡ sin α and α is a mixing angle which brings the massive neutral gauge bosons into diagonal form. Its full expression can be found in [73], but in the limit of small mixing it is well approximated by where t w ≡ tan θ w . Having the coupling (2.15) it is straight forward to calculate the decay rate into fermions, and hence find the branching ratios of the dark photon ( Fig. 6, left panel). Note we have replaced the quark masses in the charm and beauty phase space suppression factors with the lightest respective meson masses. One complication, however, is the existence of hadronic resonances between the pion threshold and ∼ 5 GeV. To gain some insight of the errors introduced, we extract the total hadronic width using the experimentally determined R Hadron factor, which is a measurement of the off-shell photon branching into muon pairs compared to hadrons [74]. Note that, unlike for the CMB constraints in [18,20], we cannot not use R Hadron directly, as we require the detailed final state photon spectrum for our calculation of the BSF limits. A comparison of the experimentally determined hadronic width to the perturbative calculation is shown in the right panel of Fig. 6. Given other uncertainties, such as the J-factors, which will enter into our limits, we deem the error introduced is acceptable for m V 2 GeV. 5

Visible photons from the decay
After having found the dark photon couplings and decay rate into SM final states we now need to find the resulting γ-ray spectrum. This is done in two steps. Firstly, for a given polarisation of V, we determine the angular distribution of decay products in the vector rest frame. We then outline how to boost this spectrum into the observer frame, where now the angular distributions of the decay products in the V rest frame converts to an energy distribution for those same final states. We now address these two issues in turn.
Example outputs from the procedure described below is shown in Fig. 7. There, we depict observer frame photon spectra for decays into electrons and b-quarks, for several parameter choices. Note the full spectrum is determined by weighting all the relevant final states by the branching fractions given in Fig. 6.

Angular distribution of V decays
Consider the angular dependence of decays of V → ff in the V rest frame. We define our coordinates such thatẑ represents the axis along which the vector is boosted in the observer frame. We are then interested in determining the distribution of decay products with respect to this axis, and accordingly define θ ∈ [0, π] to be the angle between the fermion and the boost axis in thex −ẑ plane. Taking the two circular transverse polarisations to have the explicit form µ ± = (0, 1, ±i, 0), we can determine the angular dependence as where we have defined the fermion boost (2.20) For the longitudinal polarisation, µ 0 = (0, 0, 0, 1) in the rest frame, the equivalent expression is given by In detail it is clear that the angular distribution of the fermions varies between the polarisations. When we boost to the observer frame, discussed next, this will translate into different energy distributions.

Boost of the photon spectrum
From the above, we can determine the fermion energies in the observer frame for each of the vector polarisations. However this is not the experimental quantity of interest. Instead, we aim to determine mentation in Herwig [75], specifically for the case of light dark photons. As we are using Pythia to find the final state photon spectrum, we leave the incorporation of these details relevant for mV 2 GeV for future work (see also [76]).
the distribution of photons that result from the initial hard decay V → ff . These two will coincide in the limit that the photons are produced collinearly with the fermions. Given the collinear enhancement of photon emission off a charged fermion, for certain final states this is a good approximation. For the moment let us simply assume this is true and determine the modification to the spectrum, returning to the question of when this should apply next. We define the spectrum of photon energies, E 0 , in the V rest frame as Assuming the photon is collinear with the fermions, then in the observer frame where the vector has an energy E V E D (as the initial kinetic energy is negligible) the photon energy, E, is now Importantly, we see that this energy is determined not only by the distribution of rest frame energies in Eq. (2.22), but also by the angle with respect to the rest frame, which is drawn from a distribution that depends on the polarisation of V, as determined above. In detail, and as determined in Appendix D, the spectrum in the boosted frame depends on the angular distribution p(cos θ), and takes the form where the terminals of integration are These expressions are written in terms of dimensionless quantities, in particular a boost parameter Note the absence of a factor of 2 in x arises, as after boosting in principle the photon can carry the full energy of the vector, whereas in the rest frame E 0 ≤ m V /2.

Photon spectra in the V rest frame
Equation (2.24) provides the photon spectrum in the observer frame, assuming the photons in the rest frame are collinear with the initial fermions. In this case, it is clear that the vector polarisation enters centrally through p(cos θ) (note that p[cos θ] = 1/2 corresponds to the unpolarised decays). Further, note that this result does not assume E V m V m f , as in parts of the parameter space that will not be true.
To determine the full spectra for a given set of model parameters, we will need to use this result for the appropriate set of fermions weighted by the branching fractions given in Fig. 6. Working below m V = 100 GeV, we can neglect decay to tt, however we will still need to consider six final states: e, µ, τ , q = (u + d + s)/3, c, and b. In practice we will approximate q ≈ d, as the spectra for each of the light quarks is similar. We determine the rest frame spectra for V → eeγ and V → µµγ analytically, without assuming m V m l (see Appendix E). For muons there is also a contribution  For the electron final state, we take m V = 100 MeV, and show the spectrum for transverse and longitudinally polarised V, which in this case can have a significant impact on the spectrum. For the coloured final state, we take m V = 100 GeV, and now do not distinguish between polarisations (as described in the text there is not an appreciable difference between these for hadronic final states). In both cases we show results for two dark photon boosts, γ = E V /m V . Note that for γ = 1, x = E γ /E V ≤ 0.5, and therefore in the left plot a clear transition to that regime is observed for a small boost.
from the radiative decay µ → eν e ν µ γ which is also included, following [77]. For the hadronic final state, including the τ and quarks, we use Pythia to generate the spectra. 6 With the rest frame spectra in hand, we can now revisit the question of how good an assumption it is to treat the photons as collinear with the initial fermions, as assumed in the derivation of Eq. (2.24). In particular, all of the final states above (except for the radiative decay of the muon) can be simulated in Pythia, and then boosted for each final state photon to determine the observer frame distribution. In order to simulate the distribution of initial fermion angles according to the various vector polarisations, we weight the events according to the distributions p ±,0 (cos θ) determined above.
The results of this procedure are then compared against the output of Eq. (2.24). We find very good agreement for leptonic final states, e, µ, and τ , which is unsurprising as we find the photons in this case to be predominantly collinear with the leptons. For hadronic final states, the correlation is less defined, and accordingly the collinear approximation breaks down. Nevertheless, we find that the distribution in this case is well approximated by the assumption of an unpolarised decay, p(cos θ) = 1/2 or equivalently treating the vector as a scalar.  [81,82] and Fermi [83] on DM annihilation into two photons, for bound state formation via photon emission, for two examples of E γ /m DM ratio. Also shown is the s-wave unitarity constraint for DM annihilation,

Recasting constraints on DM annihilation for BSF and level transitions
Existing constraints on DM annihilation assume that the emitted radiation has energy E ≈ m DM . Let σ ann v rel xx max @M be the maximum observationally allowed cross section for annihilation of DM with mass M into the channelX + X → xx, where X,X denote the DM particles and xx the products of the DM annihilation. If XX,XX or XX bound states form via emission of an x particle of energy E, then the corresponding constraint is found via the rescaling [58,79] where m DM is the DM mass of interest. The factor (m DM /E) 2 accounts for the different number densities of DM with mass m DM and M = E. The constraint on BSF is relaxed further by a factor 2 since only one x is emitted during BSF (in contrast to xx emitted in annihilation). Equation (3.1) applies also to exothermic level transitions that follow collisional excitations of DM bound states. In this case, σ BSF should be replaced by the cross section of the scattering process that causes the excitation, while E corresponds to the energy dissipated in the de-excitation. Note that in the case of multicomponent DM, Eq. (3.1) may have to be adjusted to account for the potentially different densities of the DM components participating in the BSF or collisional excitation processes. An example recasting for xx = γγ is shown in Fig. 8. The observational constraints come from the Planck [80], H.E.S.S. [81,82], and Fermi [83] collaborations. It is simple to repeat this exercise for different channels. As seen from Fig. 8, the constraints weaken for lower E/m DM , due to the number density factor.
In many models, however, such as the one considered in Section 2, the annihilation channel is an exotic one involving non-SM mediators that subsequently decay into SM particles. Although indirect searches have been applied to DM annihilation into exotic channels (see e.g. [17,18,20,85]), the resulting bounds are typically given in terms of a number of fundamental model parameters and are difficult to recast for the purposes of level transitions and BSF. Constraining such processes necessitates reanalysing the observational data and casting the results in terms of the energy dissipated in the transitions. In the following, we carry out such an analysis for transitions occurring via emission of dark photons decaying into SM particles.

The BSF rate and photon flux
We now return to the specifics of the model of Section 2. Outside the parameter space where dark electrons may thermalise and get ejected from the halo, we can assume that the local dark proton and dark electron densities are equal, n p = n e . Then, the total DM mass density is where f ion is the ionisation fraction (2.12). The BSF rate per unit volume is where σv rel BSF is the averaged BSF cross section (2.7). In the case of level transitions, this factor must be appropriately adjusted. Provided that the level transitions follow collisional excitation processes, then it remains true that d 2 N/(dV dt) ∝ ρ 2 DM , which ensures that the following analysis applies with the appropriate rescaling. Here we focus on BSF and shall not elaborate on the specifics of excitation and de-excitation processes.
Next we define the differential photon flux incident on the detector as where dA is an infinitesimal surface area of the detector. For a source at proper distance r only dA/(4πr 2 ) of the produced photons will reach the detector. We thus have where dN γ /dE is the visible photon spectrum resulting from BSF. Going to spherical coordinates dV = r 2 drdΩ we find where the J 0 -factor is given by and Σ is the observed area of the sky. Note in the limit m e → m p and f ion = 1 we recover, as required, the 1/16π prefactor for annihilation of non-self-conjugate DM. Going from Eq. (3.5) to (3.6) assumes that either the velocity distribution of the DM particles remains the same along the line of sight, or that (σv rel ) BSF is velocity independent. In the present case, none of these assumptions is generally true, since (σv rel ) BSF is velocity dependent as discussed in Section 2, and the DM velocity distribution within the halos varies along with ρ DM . This implies that a more refined treatment may be necessary, as we discuss next.

The J-factor velocity dependence
Let us write the BSF cross section of Eq. (2.7), as (σv rel ) BSF ≡ (σv rel ) 0 S(v rel ), where (σv rel ) 0 is velocity independent. We can rewrite the differential photon flux arising from the BSF, Eq. (3.6), to take into account the velocity dependence: where J is now the effective J-factor, which encodes the DM density, and in which the velocity dependence of the cross section has been absorbed. The full expression is [86] where f ps is the phase-space density of the dark protons and the dark electrons; since indirect signals are expected only from the regions where p and e do not thermalise, f ps is independent of the ion mass and thus the same for both species. As we have seen in a previous section dN/dE γ is a function of m V and the binding energy. 7 By using appropriate J-factors, we hope to scan over some choices of m V and the binding energy, and use Fermi-LAT data to constrain the combination of factors in the square brackets in Eq. (3.8). This factor can then be written in terms of the underlying parameters of the model and hence eventually constrain the scenario. The J-factors have been derived for S(v rel ) = v −1 rel , v 0 rel , v 2 rel , v 4 rel in Ref. [87], where the DM density and velocity dispersion were determined as functions of the radial coordinate r through a spherical Jeans analysis. Nominally these four cases are termed the Sommerfeld-enhanced (SE), swave, p-wave, and d-wave J-factors respectively. Due to the finite mediator mass, however, (σv rel ) BSF scales as v −1 rel for v rel m V /µ D , but as v 2 rel for v rel m V /µ D , as discussed in Section 2. (Note though that (σv rel ) BSF is Sommerfeld enhanced even in the latter velocity range.) To fully take this into account, we would need to re-calculate the J-factor for each choice of m V /µ D . This introduces further technical difficulties. The photon spectra depend only on m V and E V E D , which is convenient for extracting the limits on the flux, as introducing further parameters is computationally expensive. We want to avoid doing this. Furthermore the J-factors carry a large uncertainty. So we proceed by estimating the error incurred by using the pre-calculated J-factors as a simplifying approximation.
To gain some insight into this error, we can estimate the implied averaged velocity dispersion by comparing the J-factors for the different cases. If the DM density could be factored out of the velocity integral, the respective J-factors would scale as where x ≡ 2/v 2 c parametrises the velocity dispersion v c . We can extract the implied velocity dispersion, following the above assumption, by taking a ratio of J-factors. For example, the central values of the J-factor for Draco I given in [87] are  table I of [87], we find the largest discrepancy to be a factor of 0.57 for Hydrus I (SE velocity dispersion in the p-wave J-factor). We thus estimate the uncertainty introduced by neglecting the r-dependence of the velocity dispersion as a factor of a few. The reason we can do this is that if ρ DM did indeed factor out of the J-factor, i.e. there is no velocity dependence on ρ DM , then the v c would match when using the different ratios of J-factors above. By using the mis-matched velocity dispersion in the (incorrectly) factorised J-factor, we therefore obtain an estimate on the size of the effect of the r-dependence of the velocity dispersion on the J-factor. To be somewhat conservative, we will derive constraints using both sand p-wave J-factors below, which will provide further insight into the error incurred, and from which the weaker limit can be chosen.

Limits from Fermi-LAT observations of dSphs
We use Fermi-LAT observations towards dSphs to set constraints on the expected photon flux, and, ultimately, on (σv rel ) 0 in Eq. (3.8).
We use about 10 years of Fermi-LAT data, collected from 500 MeV up to 500 GeV. The data set and analysis pipeline strictly follows the procedure presented in Ref. [88]. In particular, we adopt datadriven s-wave J-factors obtained through a new dynamical analysis of dSphs which does not impose any prior knowledge (nor parameterisation) about the dSph DM density profile. A similar data-driven approach is applied for the determination of the background probability distribution function at the dSph position (we refer the interested reader to methodological details presented in [88]). To set constraints on the model under study, we use a standard profile-likelihood method by fully profiling over J-factor and background uncertainties. To improve the statistics (and sensitivity), we stack together the four most constraining dSphs (Draco, Sculptor, Ursa Minor, and Leo II), as explained in [88]. We conveniently normalise the signal using the combination f 2 ion /(f ion m p + f ion m e + [1 − f ion ]m H ) 2 = 1/(100 GeV) 2 , and we therefore set a 95% C.L. upper limit on (σv rel ) 0 . This can easily be rescaled when comparing the limit to the prediction at a given point in model parameter space.
The constraints in terms of the dark photon mass and energy are shown in Fig. 9. We provide the limits as a supplementary data file which can be used to constrain models with kinetically mixed dark photons. The constraints obtained using the s-wave J-factors apply on σv rel BSF , independently of the velocity scaling of the cross section, in the approximation where the DM velocity dispersion is nearly constant within the regions that contribute significantly to J. We also run the analysis for with prefactor normalisation f 2 ion /(f ion m p + f ion m e + [1 − f ion ]m H ) 2 = 1/(100 GeV) 2 . Note the limit on the flux becomes much stronger around m V ∼ 5 GeV due to the more efficient production of ordinary photons. The constraint on the velocity-independent part of the cross section is ∼ 8 orders of magnitude stronger when using the s-wave J-factors in comparison to the p-wave ones, as expected since v rel ∼ 10 −4 for dSphs. Note that the constraints obtained using the s-wave J-factors can be applied on the averaged σv rel , independently of the velocity dependence of the cross section, provided that the DM velocity dispersion is approximately constant within the regions of the halo that contribute significantly to the J-factors. p-wave J-factors and show the resulting limits in Fig. 9. In this case, J-factors values are taken from [87] and we model their distribution with a log-normal probability distribution function. For comparison with s-wave results, we use the same four dSphs for the stacked analysis.
These constraints apply as long as the dark photons decay within the area encompassed in the J-factors, which corresponds to 0.5 deg circle centered on the dSph galaxy under consideration. The closest of the four dSphs used in the analysis is Ursa Minor, at a distance of about 60 kpc [89]. To be conservative, we shall require that the dark photons decay within 1/10 of the corresponding radius, i.e. γcτ V 10 18 m, taking into account their boost factor at production, where g dec stands for the accessible decay channels. Note that this rough estimate neglects resonant features in the dark photon decay.

Comparison of constraints to model predictions 4.1 DM annihilation in the symmetric limit
We first use our results to constrain DM annihilation in the symmetric limit. For this we assume a standard secluded WIMP type scenario with equal number of p andp. The coupling α D is set to return the correct relic abundance [2,4]. We can then set E V = m p and include a multiplicative factor of two in the flux as each annihilation creates two dark photons and hence twice the number of visible photons as in our expressions for dN γ /dE. The limits are shown in Fig. 10. Note that on the Sommerfeld resonances, which show up as the thin constrained regions on the right of the plot, the cross section can be much larger today than at freeze-out. Shown in Fig. 11 is the limit on the cross section itself for different choices of m V .

Dark atom formation
We now apply the constraints to the atomic bound state formation in our dark sector. The constraint is given in terms of m V , and E V . The underlying model parameters are m V , α D , m p , and m e . Here we visualise the parameter space by fixing m V and E V , varying α D , m p , and choosing m e in order to return the required E V . Typical results, showing newly constrained regions of parameter space, are displayed in Fig. 12.
As can be seen, the novel constraints currently rule out only small areas of parameter space. For this reason, and considering the uncertainties on the DM velocity profile in the dSphs, we have not performed a velocity average over the DM distribution but simply set the velocity to an illustrative value from Eqs. (3.14) and (3.15), namely v rel = 20 km/s. The analysis has been performed using the limits from both the sand p-wave J-factors. With this choice of v rel the resulting constraints on the BSF cross section differ by a factor of ≈ 4. Note the condition v rel < m V /µ D is satisfied over the entire range of the plots in Fig. 12. Nevertheless, the non-trivial v rel dependence of the cross section means the assumed v rel does not entirely factor out for the p-wave constraint, as would be the case for a pure v 2 rel dependence. This shows the underlying error incurred through this approximate technique. To overcome this source of uncertainty it would be necessary to fully account for the non-trivial velocity dependence of (σv rel ) BSF when determining the J-factor from the estimate of the underlying DM phase space distribution.
Limits could also be derived using observations of the Galactic Centre, which features a higher v rel , and hence higher (σv rel ) BSF . Albeit, one must then deal with the complication of the well known excess in Fermi-LAT observations of the Galactic Centre over the standard background modelling, e.g. see [90][91][92][93][94][95][96][97][98][99].

Variations
Finally we can consider variants of the above model. For example, if there is another dark sector force in addition to the U (1) D , the binding energy of the composite state can be made larger, while keeping m p small enough to not suppress the signal due to the falling number density, and keeping α D in the perturbative range. Such a setup has recently been considered in Ref. [58]. Here, BSF occurs when the upper (N + ) and lower (N − ) components of a dark baryon isospin doublet, with opposite U (1) D charges, combine and emit a dark photon. The total binding energy is now no longer solely determined by the U (1) D but also involves an additional force, e.g. originating from a local dark SU (3) D symmetry. The cross section has been calculated in [58] and we extract it from their Fig. 3 for an example parameter point. We then confront it with our constraint from the dSphs in Fig. 13.
We also compare to the approximate dSph constraint derived in [58]. This was found by using the scaling as in Eq. (3.1). For the constraint σ N + N − →V V v rel the authors of [58] took the available limit for DM annihilating to V V followed by the 100% decay of V → τ τ from [100] and then weakened it   for three sets of parameters. The constrained regions are shown by a red contour. The grey contours show the would-be constrained regions if the limit on the flux were improved by a factor of ten. The DM relative velocity is set to v rel = 20 km/s. We have enforced m e < m p which implies we cannot reach the required E V in the white regions of the plots. The red dashed line shows the minimum allowed coupling to avoid overclosure in a standard thermal history [2,20]. The constraint using the s-wave (p-wave) J-factor is shown on the left (right). The p-wave constraint is a factor of ≈ 4 stronger. Figure 13: Constraint on the bound state formation considered in Ref. [58] from our analysis and compared to the approximate constraint derived in [58] using the results of [100]. Although the constraints here are in rough agreement, our constraints can be applied to a wider range of dark photon masses. by 1/0.1. This last factor is included as a dark photon with m V ≈ 5 GeV decays into τ τ with a branching fraction of around 0.1 (see Fig. 6). From Fig. 13 we see the constraint derived using these approximations is not too far off our constraint which takes into account the various decay channels of V more precisely. The key point is that using our results such models can be constrained more widely, with fewer assumptions, and greater ease.

Conclusions
The radiative formation of DM bound states, as well as exothermic level transitions between bound levels, provide novel sources of signals that can be probed via indirect searches. The existence of bound levels -a consequence of long-range interactions -is an important feature of many selfinteracting and/or asymmetric DM models. Unitarity arguments along with various model-dependent considerations suggest it is also a generic characteristic of (symmetric or asymmetric) thermal-relic DM in the multi-TeV mass regime and above. As our DM searches move beyond the paradigm of 100 GeV -1 TeV symmetric thermal-relic DM, identifying and exploring such novel signatures becomes essential.
In this paper, we employed indirect searches to derive constraints on the formation of DM bound states that occurs with emission of a dark photon kinetically mixed with hypercharge. We used Fermi-LAT observations of dSphs, but our analysis can of course be extended to other experiments, such as H.E.S.S., or other celestial regions of interest, such as the Galactic Centre. Our results are cast in terms of the amount of energy dissipated and the DM mass, which determines the number density of the dark particles. While the radiated energy in DM annihilation is of the order of the DM mass, BSF occurs with dissipation of a smaller amount of energy that depends on the underlying dynamics. Our results are therefore applicable to a variety of DM models where BSF occurs via dark photon emission, and reproduce also constraints on DM annihilation into dark photons. In addition, we have discussed the recasting of existing constraints on DM annihilation into SM particles, to apply on BSF.
In the course of this work we developed the treatment of a number of subtleties, namely the effects of the dark photon polarisation states, and the non-trivial velocity dependence of the bound state formation cross section. The latter means the conventionally given J-factors do not fully fit the requirements of the model. We estimated the error introduced by using an approximate technique. If limits eventually become more constraining on such models, a more careful treatment of the J-factors may become necessary. Further improvements can also be made by taking into account the effect of low lying QCD resonances on the photon spectrum produced in the cascades of the dark photon decay products [75].
We have considered and applied our constraints on a simple dark QED model of asymmetric DM that implies the existence of dark atoms forming via emission of dark photons. We determined the BSF cross section, the DM ionization fraction, and the γ-ray spectrum arising from the cascades of the dark photon decay products. The combination of these elements allowed us to predict the photon flux resulting from BSF as a function of the underlying model parameters. Thus allowing us to derive novel constraints on the parameter space. We found that the predicted flux typically lies below the derived limit, except for some resonance peaks at relatively large values of the dark coupling α D 0.1. Furthermore, variations of the model can lead to somewhat larger signal predictions [58] which we also briefly explored.
We also showed that our constraints can be applied to the annihilation and the decay products of unstable bound states of symmetric DM. In agreement with previous studies [18], we observed that in this case, the low-energy dark photon emitted in the formation of the DM bound states does not constrain the model any further due to the suppression of the DM number density by the large DM mass. However, this result does not preclude that the low-energy radiation emitted in the formation of (unstable) bound states can yield an observable signal. It is possible for example that in other DM models the spectral features of the low-energy radiation produced in BSF differ from those of the high-energy radiation emitted in DM annihilation or in the decay of unstable bound states, and render it competitive.
Crucial input in predicting the signals generated by BSF -and in fact in predicting any manifestation of DM today -is the preceding cosmological history. In the model of atomic DM considered here, the cosmological evolution determines the residual ionized component of DM that is available to form bound states today. A large BSF cross section may imply suppressed indirect signals today because the DM has already formed deeply bound atomic states in the early Universe. The details of the interplay between cosmology and phenomenology depend on the DM model and it is essential to compute these two self-consistently. For models that feature long-range interactions, the formation of stable or unstable bound states in the early Universe can critically affect all expected phenomenology of DM today [4,36,46].

A Further constraints on the dark photon
The leading constraints come from a number of sources. In Fig. 14 we have chosen to show the more stringent constraints, also including the latest supernova and BBN (Big Bang nucleosynthesis) limits, important for smaller .
• Electron g − 2. The strongest constraint in the top left corner of the plot comes from the anomalous magnetic moment of the electron [103].
Recently a stronger constraint has been set by considering energy transfer by dark photons from the centre of the supernova to the outer layers, which can affect the explosion [117].
At lower values of a constraint has been set by considering dark photons escaping Galactic supernovae [118].
• BBN/CMB. The constraints have recently been updated for particles decaying electromagnetically using a BBN code [102,119] (a comparable limit is derived from CMB N Eff measurements -late decaying dark photons do not heat the SM neutrinos, lowering N Eff [102]). Taking the limit on the lifetime from [102], we can convert it into a limit on for a given dark photon mass. Note the limit shown is actually somewhat conservative, as our dark photons will have a ∼ 40% higher temperature at BBN than what is assumed in [102] (see above). The results qualitatively match those derived analytically in [18]. Quantitatively, the results from the BBN code are somewhat more stringent. The sharp cut-off at high masses is due to the limited range of the plot in [102], although the limit is known to become progressively weaker as m V increases [18].
• Relic Dilution. Although not a constraint, we show on the plot the area in which relics are diluted by the entropy injection, due to the long lived dark photon [68]. Somewhat arbitrarily, we show a contour for which Y DM → Y DM /2 following dark photon decay.
• Direct detection. The DM will induce nuclear recoils in direct detection experiments via tchannel exchange of the dark mediator V and the SM Z boson. We derive a constraint from XENON1T data as described in Appendix B, considering only the ionised DM component. Unlike the other constraints, the direct detection limit also depends on α D and m p , which determine the scattering of dark protons on the target, as well as on m e which affects the residual DM ionisation fraction. The parameters chosen for the example shown in Fig. 14 correspond to a rather stringent exclusion contour. Analogously to the indirect detection signals, the direct detection rate does not increase monotonically with the coupling α D ; large α D may imply the efficient formation of deeply bound dark atoms in the early Universe, whose interaction with the target nuclei is partially screened due to their zero net charge.
Some overall remarks are now in order regarding Fig. 14. We conclude that even with a choice of α D and m p which results in a stringent direct detection constraint there is still unexcluded parameter space for which the dark photon does not lead to additional dilution of relic densities. Note also that since the publication of [20] the area relevant for self-interacting DM, m V O(10) MeV, has largely been ruled out from the updated BBN and SN constraints (at that time a window around ∼ 10 −10 was still open and also not excluded by direct detection). We therefore do not consider large DM self-interactions in this work. For large couplings there is also the issue of apparent unitarity violation due to the breakdown of the perturbative expansion together with the possibility of low lying Landau poles.

• Unitarity and Landau poles
The running of the dark gauge coupling is described by [120] where q is the renormalisation scale, and the β function is analogous to QED and given at two-loop level by [121][122][123][124] and n F are the number of Dirac fermions. To be concrete, for m p < q we have n F = 2, for intermediate values m e < q < m p we have n F = 1, and for q < m e we have n F = 0. The result of an evaluation of the running is shown in Fig. 15. As can be seen, we are free from low lying Landau poles provided α D (m p ) 0.2, although slightly higher values are also possible depending on the demand placed on the range for a valid EFT above m p . Note our perturbative expansion leads to apparent unitarity violation in the DM annihilation cross section for α D 0.68 [2,20]. The two constraints, no low lying Landau pole and no unitarity violation, therefore lead to roughly the same ballpark constraint on α D .

B Direct detection
Nuclear recoils are induced through t-channel exchange of the dark mediator V and the SM Z boson, as depicted in Fig. 16 (the photon has no tree level coupling to the DM). Here we do not attempt a thorough analysis of the DM direct detection, and will consider only the ionised component of DM.
Due to their neutrality, the interaction of dark atoms with target nuclei is partially screened, although it may still be significant; we refer to [40,41] for details. The spin-independent dark ion -target nucleus cross-section is given by where E R is the recoil energy, M T is the mass of the target nucleus, g V p (g Zp ) is the effective coupling of the V (Z) to the dark matter, c V p,n (c Zp,n ) is the effective coupling of the V (Z) to the SM proton and neutron, A T (Z T ) is the atomic mass (electric charge) of the target nucleus, and F Helm is the Helm form factor [125,126]. The effective couplings to the DM are given by where α is a mixing angle which brings the massive neutral gauge bosons into diagonal form whose approximate expression can be found in Eq. (2.17) (the full expression [73] is used in our code). The couplings of the vector bosons to the nucleons can be derived from Here the couplings to the chiral components of the fields are given by [73] where we remind the reader that T 3f (Y f ) is the eigenvalue of the weak isospin (weak hypercharge) of the chiral field f , and η is given in Eq. (2.16). It is well known to dark photon aficionados that in the limit m V M Z , the c V f couplings to the fermions become proportional to Q f . Furthermore, the Z exchange becomes suppressed compared to the V exchange, due to the far more massive propagator. The above cross section then reduces to an electromagnetic one suppressed by an 2 factor and modulo the finite m V mass. In this limit we may write where g EM is the electromagnetic coupling strength. This cross section has been used in a number of previous studies, e.g. [18,20]. Amusingly, the full cross section, Eq. (B.1), reduces to the same limiting behaviour also for heavier dark mediator masses. To see this, note that for the non-electromagnetic type coupling of V to be in effect, m V 10 GeV 10 MeV 2E R M T , as the recoil energy is limited by the non-relativistic velocities of the DM in the halo. We can therefore ignore the momentum exchange in the propagators. In the limit of a small mixing, the couplings can be approximated by where δ ≡ m V /M Z , and c SM Zf is the SM coupling of the Z boson to chiral fermion f , which can be found by using Eq. (B.8) and taking the appropriate limit. Substituting the above approximate forms into Eq. (B.1) and ignoring the momentum exchange, one finds the different couplings and masses associated with the two propagators simplify down to independent of δ, which is just the same as Eq. (B.9) albeit with no momentum exchange. The rate of nuclear recoils per unit of fiducial target mass is given by where ξ T is the mass fraction of the target nucleus. Here a number of astrophysical parameters enter for which we assume the standard halo model: ρ = 0.3 GeV/cm 3 is the local DM energy density, f E ( v) is the DM speed distribution in the Earth's frame, given a Maxwellian DM velocity distribution in the halo frame with peak DM speed v 0 = 220 km/s and v Earth = 232 km/s, v esc = 544 km/s is the Milky Way's escape speed, and v min is the minimum speed for which DM particles can provide a given recoil energy E R [127]. We find the limit on the model by confronting it with the latest XENON1T results [128]. Constraints from LUX [129] and PANDAX [130,131] are expected to give similar results. At low DM masses, CRESST-III [132], CDMS [133], CDEX [134], and DarkSide [135] provide more stringent constraints, see e.g. the analysis in [18,20,136]. Alternatively the Migdal effect can be exploited [137][138][139][140]. Here we shall focus on the limits for m p 50 GeV using a simple analysis. More sophisticated analyses taking into account the shape of the spectrum are of course possible [141].
Note we have multi-component DM in our model. For sufficiently heavy masses for the DM components, m p , m e 50 GeV, away from threshold effects, this increases the expected number of scattering events by a factor of two, as n e = n p . This holds provided m e m p , which we assume here, so the former component is negligible for the total DM energy density, otherwise there is a suppression as can be seen in (B.13). This factor of two is included in our limit. For a more detailed study of direct detection of multi-component DM see [142].
To set a limit we use Eqs. (B.1) and (B.13) convoluted with the best fit total efficiency of the detector, shown in Fig. 1 of [128], to find the expected number of events in XENON1T for our model given , α D , m p , and m V . The XENON1T collaboration has reported 14 events in their nuclear recoil signal reference region in 278.8 days of exposure time of their 1.3 tonnes of fiducial mass, see the second column, table I of [128]. The estimated background is 7.36 ± 0.61 events. We take the 90% C.L. limit which corresponds to DM contributing 12.8 events [143]. We find an exclusion by demanding the expected number of events at a given parameter point in our model not exceed 12.8. The result of such a procedure is shown in Fig. 14. Although this is a simplified procedure, for DM masses m p 30 GeV, it returns a limit on the generic spin-independent cross section matching that of the XENON1T analysis within a factor of two. Thus it is sufficiently accurate for our purposes here.

C Dark sector temperature
Let us denote the visible sector temperature with T and the dark sector temperature T D . Following from independent conservation of entropy in each sector, the temperature ratio is given by where τ i is the initial temperature ratio at temperature T i and h SM (h D ) count the effective entropic degrees-of-freedom in the SM (dark) sector. Prior to the decay of the dark photons, the effective entropic degrees-of-freedom in the dark sector may be modelled as where we model the disappearance of a massive species from the thermal bath with the ratio of the number density to the massless number density,ñ(x) = (x) 2 K 2 (x)/2, where K 2 (x) is the modified Bessel function of the second kind of order two. Now we wish to find τ (T ). Due to T D entering on both sides of Eq. (C.1), one can not trivially evaluate τ (T ) analytically. Nevertheless, as long as τ does not depart too far from unity, one can easily estimate it by taking into account the various mass thresholds in the dark sector, together with the SM degrees-of-freedom. To obtain a more accurate evaluation of τ (T ), an iterative approach can be used. The result of such an evaluation is shown in Fig. 17. Very similarly one can of course also find τ as a function of T D .

D Boosting spectra between frames
Here we provide a derivation of the boosted spectrum result provided in Eq. 2.24. Before converting to dimensionless parameters, the boosted spectrum can be written as where E and E 0 are the photon energy in the observer and V rest frames, respectively. In detail, we know that the energy of the photon in the observer frame is given by Eq. 2.23. This energy depends on both the energy and angle of the photon in the V rest frame, each of which are drawn from the distributions dN/dE 0 and p(z) respectively. We obtain the full spectrum by simply marginalising over both of these distributions. Converting to dimensionless quantities, we have

(D.2)
Recall B = (m V /E V ) 2 , x 0 = 2E 0 /m V , and x = E/E V . Now we will use the δ-function to perform the angular z integral. To do so, we must consider where the δ-function has support. To begin with, as x 0 ∈ [0, 1] generically, we have For the δ function to have support, we require Accordingly, we conclude which is the result quoted in the main text. Note if we are not considering polarised V decays, but just averaging over all polarisations, then we take p(z) = 1/2, and the result is equivalent to (B3)/(B4) of [144]. Similarly, in the large hierarchies limit ( B → 0), this reduces to (14) of the same work.

E Analytic results for Final State Radiation
We want to determine the spectrum of photons resulting from final state radiation of the form V → + − γ, where = e, µ. Conventionally in the literature, the form used is See, for example, (A2) of [144]. Recall here x = 2E γ /m V and l = m 2 /m 2 V . The above is an expansion in l , so it assumes l 1. Nevertheless, we are considering small vector masses, all the way to l ∼ 1, and thus this approximation will not be valid. Thus we need a more general result.
We can obtain this from the calculation in [145] for e + e − → QQg, where Q is a heavy quark. From that result, we determine From this form we can see straightforwardly, that in the limit l → 0, this reduces to Eq. (E.1) up to O( l ) corrections. The full result shown here also agrees with the calculation in [76]. From the full result, we can determine the kinematic limits on the photon energy. The minimum photon energy is 0, whereas the maximum is when the photon is emitted in the opposite direction of the + − , which are collinear and of equal energy. Then we have which rearranges to give a maximum of x = 1 − 4 l , and hence x ∈ [0, 1 − 4 l ]. We can see from the above that if x > 1 − 4 l , then 1 − 4 l /(1 − x) becomes imaginary. In the l → 0 limit, we have x ∈ [0, 1]. Numerically, we can compare the exact and approximate expressions. This is done for two different values of m V in Fig. 18 (all spectra in the V rest frame). We see that for m V ∼ 2m e there is a significant difference.