{'fraction_non_alphanumeric': 0.03222481918398145, 'fraction_numerical': 0.008999392679021662, 'mean_word_length': 4.4726558565817305, 'pattern_counts': {'":': 0, '<': 3, '<?xml version=': 0, '>': 3, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 6, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Dark Nucleosynthesis: Cross-sections and Astrophysical Signals

We investigate dark matter bound-state formation and its implication for indirect-detection experiments. We focus on the case where dark matter is a baryon of a strongly-coupled dark sector and provide generic formulae for the formation of shallow nuclear bound states on emission of photons, and W and Z gauge bosons. These processes can occur via electric and magnetic transitions, and give rise to indirect signals that are testable in monochromatic and diffuse photon measurements by Fermi and HESS. We also study the validity of factorizing the bound-state formation cross section into a short-distance nuclear part multiplied by Sommerfeld-enhancement factors. We find that the short-distance nuclear potential often violates factorization, modifying in particular the location of the peaks associated with zero-energy bound states. Finally we revisit bound-state formation of a (weakly-coupled) Minimal DM quintuplet including isospin-breaking effects, and find it gives rise to indirect-detection signals that are compatible with current bounds.


Introduction
Bound states are a critical ingredient of our universe, and of cosmology in particular. It is a natural possibility that Dark Matter (DM) particles might also bind into more complex structures. Indeed if DM is subject to attractive forces, bound states exist in the spectrum. The formation of these bound states can have consequences for the DM abundance and indirect-detection signals, as well as the formation of multi-component DM.
In our universe two types of bound states play an important role: hydrogen atoms and nuclei. While the former is weakly coupled and can be studied with ordinary perturbation theory, study of the latter requires non-perturbative methods to control the strong interactions. In this work we will study in detail the analogous problem in the dark sector: the formation of DM bound states in the strongly-coupled regime.
Many works have considered the possibility that DM is charged under a new abelian gauge symmetry, which allows the formation of hydrogen-like DM bound states, see [1] for a review. In [2] we undertook a systematic study of hydrogen-like bound states that arise in any weakly-coupled abelian or non-abelian gauge theory (see also [3]). This effect is relevant for WIMP DM candidates such as Standard Model (SM) electroweak multiplets, as it can significantly modify the annihilation cross-section that determines the DM thermal abundance. For example we found that bound state formation increases the mass at which an SU(2) quintuplet reproduces the critical abundance from 9.5 TeV [4] to M ≈ 14 TeV. It also leads to novel indirect-detection signals due to the emission of quanta in the formation process.
The production of strongly-coupled DM bound states was first studied quantitively in [5]. In this work residual strong interactions support the formation of more complex structures. The canonical example is a scenario where DM is a baryon of a confining dark sector [6] (see [7] for a review). The dynamics is similar to that of QCD, so nucleons and nuclei may be likewise stable due to a conserved (dark) baryon number. Computing the relevant nuclear cross sections seems unfeasible due to the strongly-coupled nature of the problem. Fortunately, as was understood and exploited in nuclear physics, this is not the case if the bound states are shallow. As shown by Bethe and Longmire in the '50s [8] and systematically rederived using nuclear effective theories in the '90s [9,10], the phenomenology of the bound state formation depends only weakly on the details of the underlying potential once the binding energy is fixed. This observation allows us to reliably compute the cross section in a controlled effective-theory expansion. These results can be applied to dark sectors with strong interactions such as [11][12][13][14][15] In this work we refine and extend previous studies in several directions. First we address the issue of factorization of the cross section into long-distance effects (Sommerfeld enhancement) and a short-distance nuclear cross section. We show that standard factorization is often violated because the zero-energy bound states associated to Sommerferld peaks are modified by the short-distance potential. This effect can be taken into account in the cross-section computation by including both long-and shortdistance contributions to the potential in computing the wavefunctions that determine the overlap. We present general formulae for bound-state formation by electric and magnetic emission of electroweak gauge bosons.
Next we consider indirect-detection signals associated to bound-state formation [16] (see [17][18][19] for weakly-coupled realizations). If DM is a dark baryon, dark deuterium can be formed through emission of quanta of energy equal to the binding energy of the state. In the simplest case where DM is also charged under the SM this is automatically realized through the emission of electroweak gauge bosons. This process gives rise to monochromatic photon lines that are constrained by observation of the Galactic center, as well as a diffuse photon signal from emission of W and Z gauge bosons and their subsequent radiation. For example, when the DM baryon is a triplet of SU(2) L the cross section can be within the reach of current experiments such as FERMI and HESS in the vicinity of its peaks. We also discuss bound-state formation of Minimal Dark Matter [20], and find that the associated signals are consistent with the current bound on the thermal mass from relic density considerations. This paper is organized as follows. In Section 2 we present our main results and illustrate them in the context of a simple U(1) toy model. In Section 3 we provide general formulae for the formation of two-particle bound states by emission of a gauge boson. In Section 4 we apply our formalism to the simplest dark baryon model: a complex SU(2) L triplet. We first compute the rate of cosmological production of 'dark deuterium', improving previous estimates by including isospin-breaking effects, and then derive astrophysical signals due to bound-state formation. In Section 5 we apply the same tools to Minimal DM: an elementary quintuplet of SU(2) L . We summarise the result in Section 6. A series of technical appendices follow.

DM DM
DM bound state gauge boson Figure 1. DM bound-state formation. When DM is charged under a long-range force the initial state wavefunction is strongly distorted before the short distance nuclear process, depicted as red blob, takes place.

Formation of shallow bound states
We are interested in computing the formation of bound states, which we will refer to generically as 'dark deuterium', via emission of SM gauge bosons. This is always kinematically allowed for photons, while for SU(2) gauge bosons the binding energy must be larger than their mass in order for the process to take place at rest. For the purposes of this work we take DM to be a spin-1/2 dark baryon triplet, V , of SU(2) L . This could arise, for example, from an SU (3) D gauge theory with fermions in the fundamental representation of dark color, which are also triplets of SU(2) L [6]. The lightest nucleon is the neutral component of the baryon triplet, which has a 165-MeV mass splitting from the charged component due to EW-loop effects, and it is stable due to dark-baryon-number conservation. Since DM has no electric charge, bound-state formation by emission of photons cannot take place at tree level, and requires the inclusion of longdistance effects. If DM is an electroweak multiplet, Sommerfeld enhancement (SE) due to multiple exchanges of electroweak gauge bosons can turn a pair of V 0 V 0 into V + V − . The latter can then emit photons on bound-state formation. Alternatively the strong interaction can also allow for the same transition.
The problem of computing bound-state formation in the presence of long-distance effects such as SM gauge interactions is complicated by the presence of multiple scales. We wish here to clarify the issue of factorization of the total cross section, which will often not hold for bound state formation 1 .
We start by listing the relevant scales of the problem: • r dB = 1/(M v rel ): the de Broglie wave-length is typically the largest length scale in the problem where the particles can be treated as free.
• r ∆ = 1/ √ 2M ∆: the length scale where the splitting between charged and neutral components (∆ =165 MeV) becomes important. For energies below 2∆ the V + V − pair does not exist as an asymptotic state, and its wavefunction goes to zero at distances larger than r ∆ . Below this distance scale the the system is, to a good approximation, isospin-symmetric, and the wavefunction of the charged pair is related to that of the neutral pair by a Clebsch-Gordan coefficient.
• r W = 1/M W : SU (2) interactions become relevant and approximately symmetric at distances smaller than r W . This the scale where the SE starts to build up.
The 'Bohr radius' is the distance at which the SE approaches the asymptotic value commonly used in the literature.
• a i,f : the scattering lengths of the initial and final states of the strong nuclear interaction. This corresponds to the size of the shallow bound states in the initial and final channel.
• r π = 1/M π : the range of the strong interactions. For shallow bound states such as nuclei, a i,f > r π and detailed knowledge of the microscopic physics becomes unimportant. Indeed, to leading order the cross section only depends on the scattering length, as first pointed out by Bethe [8].
To first approximation the DM annihilation cross section is assumed to factorize into a short-distance part and a long-distance part as follows: where SE is the Sommerfeld enhancement encoding the distortion of the initial wavefunction due to exchange of light mediators, while σ short is computed from microscopic physics. This factorization is justified in the case of heavy DM, because the annihilation effectively takes place at a distance r ann ∼ 1/M , the wavelength of the final states, which is much shorter that r dB [22,23]. For boundstate formation there are two (related) effects that may invalidate the factorization. First, for some regions of the strong-sector parameter space, the bound state cannot be treated as point-like. Second, for massive mediators the cross section has peaks that are associated to zero-energy bound states supported by the long-distance potential. These bound states at threshold are very fragile and can be modified by the short-distance potential. Although the latter would also spoil factorization in the case of DM annihilation, the effect is unavoidable in bound-state formation if the nucleons have SM charges. We describe each of these effects in greater detail below.
Finite size effects: The failure to factorize due to finite size effects is already evident in the production of hydrogen. For the ground state, in the electric dipole approximation and at low velocities one finds (see [24]) This formula is written as a p-wave SE times a short distance cross section that depends only on the final state. The factor e −4 ≈ 0.018, arising from the non-trivial overlap of the wavefunctions of initial and final states, encodes the deviation from the factorized limit. More generally, if the initial state feels an attraction of strength λ i α while the bound state is associated to a coupling λ f α the deviation from factorization is given by e −4λ i /λ f , so that in the limit λ f λ i (small Bohr radius) the factorization is recovered.
Similarly, for the formation of deuterium at low energies one finds [8] (σv) pn→D+γ = κ 2 where a i ≈ −23 fm is the scattering length of the initial 1 S 0 channel and γ f = √ E B M ≈ (5fm) −1 is the inverse scattering length of the 3 S 1 final state channel. In the formula above the first term in the parenthesis corresponds to the tree-level process while the second term encodes long-distance nuclear effects associated to the large scattering length of the initial state, see [25].
It is clear that factorization should be a good approximation if the size of the bound state is smaller than the range of the long distance interactions. This corresponds to, When this condition is violated the effect of long-distance physics cannot be captured simply by the value of the wavefunction at the origin and a full quantum-mechanical computation is required.
Location of the Sommerfeld peaks: The other effect that violates factorization relies on the presence of a short-range attractive potential in addition to a Yukawa-type potential for the initial state. As we will see this effect is very generic. It has long been known that zero-energy bound states supported by a Yukawa potential in the initial state give rise to peaks in the cross section at low velocities. The positions of these peaks correspond to a set of critical masses M * determined by the equation where E W is the binding energy of the bound state supported by the Yukawa potential (induced for example by electroweak interactions). This tuned condition is potentially sensitive to deformations of the potential due to other effects. Indeed as we will see, any short-range potential generically contributes to the binding energy of shallow bound states such that the zero-energy bound state appears at a different value of the mass. This violates the naive factorization of long distance and short distance effects. One way to take this into account would be to compute the SE using the full potential, including both long-and short-distance contributions. While this computation reproduces the positions of the peaks, it does not take into account the fact that the effect of the short-distance potential is already partially included in the hard scattering cross section. It is possible to improve this estimate by evaluating the SE not at the origin, but at a distance of the order of the scattering length. Clearly the safest way to proceed is to compute the cross section directly using explicit wavefunctions and evaluating the overlap between initial and final state. Importantly, while this procedure requires the choice of an explicit potential, in the limit of shallow bound states the result becomes independent of the details of the short-range potential.

U(1) Toy model
In order to illustrate these effects we first consider a simple U(1) model. In this scenario DM is composite i.e subject to a short-range nuclear potential, and coupled to a vector boson of mass M V . The static potential between two DM particles with opposite charges can be approximated as where α D is the strength of the dark U(1), and r 0 < 1/M V is the range of the nuclear interaction. As discussed in [16], strongly-coupled models with a dark photon interaction can be effectively parametrized with this type of potentials. When the nuclear potential is switched off (V N = 0) the Yukawa potential alone supports a zero-energy s-wave bound state for Let us now turn on an attractive short-distance nuclear potential (V N < 0). As we increase |V N | the zero energy bound state becomes deeper until For such a value the nuclear potential is not a small perturbation on the original Yukawa potential, but rather the converse. Indeed when this condition is met the nuclear potential alone supports one bound state. The long distance potential gives a shift to the binding energy of the nuclear bound state that can be estimated as where we used the fact that the size of the nuclear wavefunction is comparable to the scattering length, which in turn is determined by a −1 ≈ √ M E B . This shows that corrections to the nuclear binding energies are small, possibly exponentially suppressed.
We can then ask about the fate of the zero-energy bound state supported by the long distance potential, and perturbed by the nuclear potential. Naively one might think that the attractive shortdistance potential should make the zero-energy bound state more bound. However one must keep in mind that this bound state would be the first excited level of the full potential. On general grounds the wavefunction of the ground state has zero nodes while that of the first excited level has one. This implies an extra positive contribution to the energy from the kinetic term. Indeed in our example the bound state is less bound, shifting the peak of the SE to larger masses.
In Fig. 2 we plot the binding energies of the first two bound states of the full potential (6), as well as those of the pure Yukawa potential. We see that in the presence of a nuclear potential the zero energy bound state is shifted to larger masses. This implies that the Sommerfeld peaks of the cross section will appear at different values of the mass.
For p-wave zero modes the situation is different, the condition for zero-energy p-wave bound states in a spherical well is V N M = π 2 /r 2 0 ; hence for a shallow s-wave nuclear bound state there is no p-wave nuclear bound state. A zero mode p-wave state arises solely due to the Yukawa interaction, and it has the lowest energy level in this sector. Equivalently, since the wavefunction vanishes at the origin the effect of the short-distance potential is small, and the location of the peaks will not change much. This intuition is confirmed by numerical computations Let us now discuss the production of nuclear bound states. For simplicity we focus on magnetic interactions which allow the formation of an s-wave bound state from an s-wave initial state by emission of a dark photon. As reviewed in Appendix C a shallow s-wave bound state gives rise to a SE factor, where a −1 0 ≈ √ M E B is the s-wave scattering length, and p = √ M E is the momentum of the incoming state (for a generic potential V N would be replaced by the typical depth of the potential). These peaks are visible in the total cross section for bound-state formation in the U(1) toy model. We compute the cross section using the full quantum-mechanical overlap, and using the factorized cross section. Our results are shown in Fig. 3. The true cross section is given by the solid blue curve. The factorized result has a similar shape, with a normalization and shift that depends on how the Sommerfeld factor was computed. Including only the long-range Yukawa potential in the SE computation results in a shift to lower masses (orange dotted curve) as a consequence of shifted zero-energy bound states. The correct location of peaks is captured by computing the SE with the full potential, containing both short-range (nuclear) and long-range (Yukawa) components (6). A very good approximation to the true cross section is obtained by evaluating the SE at r = a 0 (blue dashed) rather than at the origin (blue dotted). Note that we only show explicitly here the breakdown of factorization due to a shift in the zero-energy bound states. In Fig. 3 we show the equivalent effect for a dark SU(2) L triplet model, which we consider in greater detail in Section 4 below. In this case the presence of several channels also leads to a dip where the cross section is strongly suppressed.  The DM potential consists of short-range (nuclear) and long-range (Yukawa) components, see (6). The binding energies with the nuclear potential turned off are shown for comparison (dashed lines). Right panel: Bound-state formation cross section (solid line) for composite U(1) model, showing shift of zero-energy bound states due to the short-range nuclear potential. The factorized cross section is shown for comparison, with the SE computed for (i) a long-range Yukawa potential only (orange), (ii) Yukawa + shortrange nuclear potential, evaluated at the origin (blue dotted) and (iii) Yukawa + short-range nuclear potential evaluated at r = a 0 , the size of the initial bound state (blue dashed).

Cross sections for bound-state formation
In this section we provide generic expressions for bound-state formation, including explicit formulae for shallow bound states. This extends the work in [2], where the formation of hydrogen-like bound states through electric dipole transitions was studied for general weakly-coupled non-abelian gauge theories.
As we discussed in the previous section the standard approach to computing the relevant cross sections in the case of bound state annihilation, i.e. convolution of the analytic short distance cross section with the SE due to long-range forces only, does not work generically. This is due to the shift of the spectrum of zero-energy bound states, as well as finite-size effects in some region of parameters. We will thus take a numerical approach where we compute the wavefunction overlap using the full potential. While this requires a choice of nuclear potential, the results are only weakly dependent on this choice in the regime of shallow bound states that will be of interest to us here. For simplicity we will parametrize the nuclear potential with a spherical well.
A note on calculability is in order. The possibility to compute an intrinsically strongly-coupled process such as nuclei formation relies on the smallness of nuclear binding energies, allowing us to describe bound state formation within a local effective theory where the pions have been integrated out. The expansion parameter is In the SM the expansion parameter is approximately 0.3, hence higher order corrections are of the order 10%, as seen experimentally.
It remains an open question whether the mild tuning of the binding energy is accidental or a robust feature of nuclear interactions. The latter possibility is supported by lattice studies [26].

Quantum mechanical computation
To compute the cross section in quantum mechanics we follow [2], to which we refer the reader for further details. The general process we wish to consider is bound state formation through emission of a spin-1 field V a with mass M a and momentum k: Here the indices i, j, i , j run over the various components of DM in the initial and final states, B i j is the bound state, and a runs over the generators of the adjoint. In the center-of-mass frame the energy of the emitted quantum is As usual we can decompose the initial-and final-state wave function in states of fixed orbital angular momentum where Y m are spherical harmonics. The radial wave function u (r) satisfies where E = M β 2 . In general this is a matrix-valued equation in the two-particle Hilbert space |i ⊗ |j and must be solved numerically. Details of our method, including the boundary conditions and normalizations are given in Appendix B.
We will include in the potential a Yukawa-like long-range component induced by SM interactions (c.f. the ubiquitous wino in the MSSM) as well as a short-range nuclear potential, as in Eq. (6). For the latter we will use a spherical well with parameters that reproduce the binding energies we select for the nuclear bound states.

Magnetic transitions.
For strongly-coupled bound states, such as the proton and neutron in the SM, the largest cross section at low velocity is usually due to magnetic dipole transitions. This observation also holds true in the context of dark sector bound states, and it can have important cosmological consequences [5]. The magnetic emission is described in general by the effective Hamiltonian in the two-particle space, where the magnetic field B ka = klm F a lm and σ k denotes the Pauli matrices. The matrices T and T are generators in the representations of the initial DM states 1 and 2 respectively, of mass M . g is the gauge coupling, and the coefficient κ is O(1) for strongly-coupled bound states.
For transitions between s-wave states by emission of a gauge boson V a with mass M a one finds, with the group theory factor, For shallow bound states, in the absence of any long-range interaction, the integral above can be computed using the effective range expansion [8,27] and is independent of the details of the nuclear potential to leading order. As shown in Appendix C the resultant short distance cross section in the low velocity regime is, where a r and 1/γ r are the scattering lengths associated to the initial-and final-state channels. 2 The same result can be derived using nuclear effective field theory techniques, see [25]. 3 Electric transitions.
If the bound states are minimally coupled to the mediator an electric dipole interaction is unavoidable. Computation of the cross section proceeds as in the magnetic case and the analogous formula for the formation of an s-wave bound state through an electric dipole interaction reads [2], (22) where the second term in the integral is associated to non-abelian interactions.
As for magnetic transitions explicit formulae can be written in the symmetric limit for shallow bound states, in the absence of long-distance interactions. Using the results in Appendix C one finds, where we neglect the non-abelian contribution as it is typically suppressed for strongly-coupled bound states. Note that in the electric case the cross section is not enhanced by a large scattering length in the initial state.
Equations (17) and (22) are completely general and can be used to compute the bound-state formation cross section for arbitrary potential and gauge groups. With these tools at hand we proceed to study the cosmology and astrophysics of DM bound-state formation, within the context of two examples that highlight the relevant physics processes in the strongly-and weakly-coupled cases.
4 Strongly-coupled model: Weak-triplet dark baryon The simplest model of nuclear DM is an SU (3) D gauge theory with fermions in the triplet representation of SU(2) L [6]; DM is the neutral component of an isospin-triplet (dark) baryon V . If it is a thermal relic its abundance is determined by its pair-annihilation into dark pions. Estimates of the rate for this process, substantiated by real QCD data, indicate a thermal mass around 100 TeV. For smaller where in the second step we took the non-relativistic limit. This interaction has the same structure as magnetic dipoles so the same selection rules apply: ∆S = 1, ∆L = 0, ∆I = 0. Hence the nuclear cross section for pion emission is naively similar to that for a magnetic transition (20), but with σ0 → 1/f 2 . For pion emission however √ M EB > Mπ, so the effective theory of nucleons cannot be used and a full strongly-coupled computation is required. If allowed this process is expected to give the largest production cross section. masses the correct relic abundance can be realized instead through an asymmetry in the dark sector completely analogous to that in the visible sector, with the DM made up of the lightest neutral baryon. We will assume here that DM is asymmetric, in which case the mass can be taken as a free parameter. We emphasize however that the signals discussed here exist even in symmetric scenarios.
The dark nuclear forces will likely give rise to larger nuclei too, with the lightest nucleus in each baryonic sector being absolutely stable if M nucleus < nM nucleon for a nucleus composed of n constituent nucleons, due to dark baryon number conservation. We will focus here on production of the nucleus with baryon number 2, or 'dark deuterium', D. Its isotopes can be decomposed by their weak-isospin representations as, Like for the wino, electroweak loop corrections will give rise to a small mass splitting between the charged and neutral components of the dark baryon V , ∆ ∼ 165 MeV for M ∼ 3 TeV. Although this splitting plays no part in the DM cosmology at early times, it becomes important at late times as the density of the charged component becomes negligible.

Cosmological abundance
The cosmological abundance of dark deuterium was estimated in [5], using the assumption of factorization for the Sommerfeld-enhanced cross section. Here we improve on our previous computation by evaluating the cross sections numerically, including the full nuclear and electroweak potentials in our determination of the initial-state wavefunction.
We assume for simplicity that DM is asymmetric with the dark baryon V as the sole component, and that only the singlet and triplet isotopes of dark deuterium D 1,3 are bound. The DM abundance can be computed by solving the following Boltzmann equation for the dark-deuteron formation process V + V → D + X, where X stands for an electroweak gauge boson γ/Z/W in equilibrium with the SM thermal bath,ṅ where n V and n D are the densities of dark baryons and dark deuterium, respectively, and (σv) eff is the inclusive thermal cross section for dark deuterium production. In the presence of several isotopes the equation above describes the total deuterium abundance with the identifications  For each value of the mass we compute the cross section numerically as a function of velocity and then integrate the Boltzmann equation above, see [5] for more details. The dark deuterium production cross section due to electric and magnetic transitions, as well as the dark deuterium mass fraction for various choices of parameters, are shown in Fig. 4. A significant fraction of dark deuterium is produced only for binding energies larger than in ordinary nuclear physics. For small binding energies the abundance is, moreover, smaller than estimated in [5].

Indirect Detection
Since DM today consists of the neutral component of the dark baryon triplet, only the two-particle DM subsector with zero charge will be relevant for indirect detection. The radial wavefunction is described by a two-component vector, corresponding to the V + V − and V 0 V 0 two-particle states, respectively. The s-wave initial state has spin 0 while p-wave state has spin 1. In both cases the static potential has the form, where A = (α em +α 2 c 2 W e −M Z r )/r, B = α 2 e −M W r /r with c W = cos θ W are the electroweak contributions, and ∆ = 166 MeV is charged-neutral mass splitting. V 1,5 N are the nuclear potential in the singlet and quintuplet isospin channel. In our numerical study we approximate the nuclear potentials with a spherical well with a depth chosen to support bound states of given binding energies, see Appendix C.
Since electric and magnetic interactions both violate isospin by one unit, isospin selection rules allow us to form an isospin-triplet nuclear bound state (spin 1) from an isospin singlet or quintuplet in the initial state. The electrically-neutral component is supported by the potential, The triplet bound state is approximately isospin symmetric for the choice of parameters we consider. In isospin space the wave function is determined by Clebsch-Gordan coefficients |10 = 1 the radial wavefunction R B (r) can be conveniently chosen to be real (see Appendix C for full form of bound-state wavefunction).
In the basis of Eq. (17) the initial-and final-state wavefunctions have the form: where R ( ) for = s, p are the radial wavefunctions for the initial-state s-and p-wave states, while R B (r) is the radial wavefunction of the isospin 1 bound state.

Photon lines
Let us now discuss astrophysical signals from photon emission on dark deuterium formation, and their measurement, either in the galactic center or in more distant dwarf galaxies. Since DM today is highly nonrelativistic (β = 10 −3 and β = 10 −4 , respectively) the photons emitted are monochromatic up to energy resolution effects, with E γ = E B , the binding energy of the final state. Recall from Table 1 and the subsequent discussion, that selection rules allow for the formation of the neutral isospin-triplet spin-1 dark deuteron D 0 3 . This state can be formed by a magnetic transition from an s-wave initial state, or equivalently by an electric transition from a p-wave spin-1 initial state. From eqs. (17,22) (T a = −T a = J 3 and g N = 2) the cross sections for these processes are: where we neglected the electroweak non-abelian contribution in the electric transition as it is suppressed for small nuclear bound states. In Fig. 5 we present the cross section for bound state formation by monochromatic photon emission from the galactic centre (for bounds from dwarf spheroidal galaxies see [30]). Unlike photons produced from DM annihilation, the energy E B of these photons is independent of the DM mass M and so all experimental bounds for bound-state formation must be rescaled according to This rescaling takes into account the reduced DM numerical density compared to annihilation into gauge bosons of DM with M DM = E B and that a single gauge boson is emitted in bound state formation. We see in the plots some generic features discussed in Section 2. The peaks of the cross section curves are related to those in Sommerfeld-enhanced DM annihilation processes, but for the magnetic transition are shifted to higher masses. The magnetic channel also features dips. These can be understood as the effect of negative interference. As shown in Appendix A the short-distance nuclear annihilation matrix has non-zero off-diagonal components, and electroweak SE allows V 0 V 0 → D 0 3 γ and V + + V − → D 0 3 γ to interfere destructively. Indeed the analogous effect is also realized for annihilating DM when the reactions proceed through several channels, see [4,31].  Figure 5. Cross section for monochromatic photon production for the process V 0 V 0 → D 0 3 + γ via the electric and magnetic transitions for different choices of nuclear potentials. The gray (red) region is the exclusion due to γ-ray lines from FERMI [28] (HESS [29]) from our galactic centre, assuming an NFW DM profile. We take β = 10 −3 .
Finally, as mentioned in [16] D 0 3 is likely not the lightest dark-deuterium isotope. As a consequence a second monochromatic photon could be emitted in transition to the ground state D 0 1 , with energy E γ = E B 1 −E B 3 . If this was the only decay channel kinematically allowed to the dark deuterium triplet state, its rate would be identical to that of the primary photon emission, giving rise to a smoking-gun signature of bound state formation.
Current experiments are already sensitive to formation of dark deuterium in the neighborhood of the peak of the of the cross section, assuming an NFW DM profile. Future experiments such as HERD [32] will improve the bounds and yield a better energy resolution for the lines.

Diffuse photons from W/Z emission
For binding energies larger than M W,Z bound-state production could also proceed through emission of W and Z gauge bosons. The cross section for these processes can be simply obtained from Eq. (17), (22), with the use of the relevant couplings and generators (T a = −T a = J 3 for Z emission, T a = −T a = J + for W ). One finds, The inclusive W cross section is twice the value above.

Magnetic
Electric Figure 6. Cross sections for W and Z emission via formation of a dark-deuterium bound state, in electric and magnetic channels. The resultant γ ray spectrum is constrained by diffuse γ−ray searches in dwarf spheroidal galaxies by FERMI [33].

Gauge boson production cross section
In Fig. 6 we compare the cross section for W and Z emission on production of a dark-deuterium bound state, with the experimental bound on diffuse photon emission from dwarf spheroidal galaxies by FERMI [33]. Close to the peaks the magnetic cross section can give rise to a signal that within the sensitivity of the experiment; this behaviour is consistent with expectations from [34].

Weakly-coupled model: Weak-quintuplet Minimal DM
The tools developed in this paper can also be used to study signals of bound-state formation in weaklycoupled models. In this section we apply these techniques to minimal DM, in the form of a fermionic SU(2) L quintuplet with zero hypercharge [20]. Including the effects of bound-state formation in the computation of the minimal DM relic density gave rise to a significant increase in thermal relic mass, from 9.5 TeV [4] to 14±1 TeV [2]. We expect existing indirect-detection bounds on minimal DM to be similarly modified on correctly accounting for bound-state formation, both in the annihilation channel, which results in emission of gauge bosons with energies equal to the DM mass, and the bound-state formation channel, where the emitted gauge bosons have energy equal to the binding energy of the final state. We update the existing phenomenological constraints in this section.

Annihilation
The annihilation rate for the neutral component of the quintuplet, χ 0 , is strongly modified by inclusion of SE in the initial state. Although a detailed computation can be found in [36], we extend their analysis to the larger range of DM masses that are thought consistent with updated relic density considerations. As in the previous example, since DM is a Majorana fermion the relevant sector has spin-0 and is electrically neutral, (Q = 0, S = 0). It encapsulates mixing between the following two- The experimental bound are extracted from Galactic center HESS data assuming an NFW profile [35]. Right panel: Annihilation into WW leading to a diffuse photons signal. We use here the same bounds as in [36] from the galactic center.
particle states: {χ ++ χ −− , χ + χ − , χ 0 χ 0 }, and the associated electroweak potential in this basis is We compute the total cross section using factorization, and including full isospin-breaking effects in the SE as detailed in Appendix A. The cross sections for annihilation of DM into pairs of (monochromatic) photons of energy M , as well as to W bosons are displayed in Fig. 7. Superposed on these plots are the bounds on monochromatic photons from the galactic center from HESS, as extracted from the 10-year line search [35], and the NFW limit from diffuse photons searches by Fermi, as extracted from Fig. 4 of Ref. [36]. We see that for M = 14 TeV where the thermal abundance of DM is nominally reproduced the signal is bordeline but still consistent with γ−ray constraints even with the aggressive NFW profile. Different values of the mass might be however excluded.

Bound-state formation
In this section we consider the indirect signals of a minimal DM quintuplet pair undergoing boundstate formation. A preliminary study in the SU(2) L -symmetric approximation can be found in [2]. Our approach allows us to compute the relevant rates without relying on the approximate symmetric limit, but instead taking into account the full potential including electroweak splittings in the initial state Bound state formation in this case has two potential paths to an indirect-detection signature. The first is due to the SM gauge boson that is emitted as a by-product, and the second is directly from the decays of the unstable bound state produced. However, since the main production channel is to a spin-1 bound state that subsequently annihilates to the SM Higgs and fermions, no significant bounds are obtained from this channel.
The leading transition mechanism for weakly-coupled bound states is via the electric coupling to electroweak gauge bosons. Recall that electric dipole interactions imply the selection rules ∆S = 0, ∆L = 1. Hence s-wave bound states are formed from a p-wave initial state while p-wave bound state are formed from s-wave and d-wave initial states. Decomposing the DM initial state in eigenstates of isospin Given that in the dipole approximation ∆I = 1 it follows that only the triplet bound state can be produced (the septuplet is not bound). The bound states produced in the process are approximately coulombian and can be easily computed and classified by their isospin I, spin S, and angular momentum (n, ) [2]. For the critical mass M * = 14 TeV the binding energies of the relevant states are given by Name I S n E B (M * )/GeV Produced from 1s 3 3 1 1 0 100 In the notation of eq. (22) the initial wavefunction reads, where R

Photon lines
From the above decomposition we can now compute bound-state formation with the emission of a photon. The energy of the emitted photon is equal to the binding energy of the bound state, which is fixed by the DM mass and the choice of quantum numbers in the initial and final states. For thermal relic mass for instance the binding energies, as read from (37) are 25 and 100 GeV. The cross section for the formation of an s-wave isospin-triplet (spin-1) bound state, and to use it we need to solve the p-wave Schroedinger equation with the potential given in Eq. (34). For completeness the analogous cross section for the formation of a p-wave isospin-triplet (spin-0) bound state is The monochromatic photon cross section due to bound-state formation in these two channels are shown in Fig. 8. For a quintuplet mass of 14 TeV, bound state production does not yield a significant constraint even in the larger p → s channel. We note however that the cross section is very sensitive to the mass motivating a precision study.

Diffuse photons
The binding energy of the s-wave isospin-triplet bound state is sufficiently large to also emit W and Z bosons in a sizeable interval around the mass where the thermal abundance is reproduced. The cross section for the latter can be obtained by rescaling Eq. (40) by the coupling strength and kinematic emission factor for the Z boson: For W emission we have instead: (43) The cross sections are reported in Fig.(9). For the nominal value of the thermal mass the bounds do not give significant constraints.

Conclusions
In this work we studied the formation of DM bound states and their potential impact on astrophysical signals of DM. Bound state formation is expected to be generic in theories where DM has strong interactions, for instance when DM is the lightest baryon of a confining dark QCD-like gauge theory. Our results are complementary and generalize in the strongly coupled regime Ref. [2] where production of perturbative hydrogen-like bound states was studied in weakly-coupled non-abelian gauge theories.
The physics of strongly-coupled DM bound states is analogous to that of deuterium in the SM, where the smallness of the binding energy allows us to compute the relevant cross sections using the effective range expansion, without a detailed knowledge of the underlying nuclear potential. If DM is charged under the SM electroweak group, bound states of DM can form through the emission of a monochromatic photon, or electroweak gauge bosons if kinematically allowed.
At the technical level we showed that the presence of a short-range (nuclear) component to the potential often invalidates the ubiquitous factorization of the cross section into a Sommerfeld factor, which encodes the long-range effect due to the electroweak potential, times a short-distance cross section. This failure obviously arises due to the finite extent of the bound states, but more importantly is due to the modification of the spectrum of zero energy bound states by the short distance potential. The latter effect would also impact DM annihilation if the DM has a strong short-distance potential. In light of this the relevant cross section for bound-state DM should be safely computed using explicit wavefunctions in order to reliably capture the interplay between the 'short-distance' and 'long-distance' physics.
After providing general formulae for bound-state formation, we considered in detail the minimal model where DM is the neutral component of an SU(2) L -triplet baryon. Formation of dark deuterium gives rise to monochromatic and diffuse photon signals can be constrained by existing FERMI and HESS measurements, and will be further tested by future experiments such as HERD.
With the tools developed in this work we also revisited SU(2) L -quintuplet Minimal DM, updating previous studies. Besides its direct annihilation to SM particles, bound state formation leads to novel signatures for this DM candidate that are unavoidable and complementary although currently consistent existing measurements for M=14 TeV.
Our formalism, which allows us to compute the formation of shallow bound states on the emission of light quanta in the strongly-coupled regime, is easily extended to other scenarios. For example if Yukawa interactions with the SM Higgs exist, bound-state formation can proceed through Higgs emission. Emission of light particles coupled to the SM at sufficiently low velocities could also be studied in a similar way. We leave these and other questions to future work.

A Factorized cross section
In this appendix we provide explicit formulae for formation of shallow bound states that are applicable when the long distance effects can be factorized in the Sommerfeld enhancement factor. In this case the cross section can be cast in form similar to the annihilation of DM. By solving the Schroedinger equation for the initial state without strong interactions we can define the matrix, where ψ 0 a is the initial state free wavefunction in the channel a. The matrix A ia encodes the Sommerfeld enhancement and is identical to one computed for annihilation.
The full cross section for the process a → i → f reads [22] where Γ f is a generalized the cross section for the nuclear process and c a = 2(1) for Majorana (Dirac) particles. Note that diagonal components of Γ f describe physical cross sections while the off-diagonal components corresponds to interference effects.
Since the initial states are charge eigenstates it is useful to give the annihilation matrices in that basis. For magnetic transition the generalized cross section matrix appearing in (45) can be written as [4], The group theory factors take into account that charge eigenstates are not mass eigenstates from the point of view of strong interactions.
For electric cross sections instead,

A.1 Weak triplet
As an example we here report the generalized cross sections matrices for the production of the s-wave nuclear bound state D 3 from an initial state with two neutral nucleons V 0 in the triplet representation of SU (2) L . For the magnetic transition the initial state is s-wave with spin-0. One finds,

B Sommerfeld enhancement and phase shifts
In this appendix we review the computation of Sommerfeld factors. Since we only consider rotationally invariant interactions, the angular and radial part of the wavefunction separate from each other. The scattering process, at large distances, is characterised by the asymptotic expression where we fix the axis of the incoming plane wave to be p · r = pr cos θ = pz. In general the plane wave along z can be decomposed as e ikz = i Y ,0 (θ) 4π(2 + 1)j (kr), where j are the spherical Bessel functions and Y ,m=0 (θ) = (2 + 1)/4πP (cos θ) the spherical harmonics. Therefore, by exploiting the behavior j (ρ) → sin(ρ − π/2)/ρ for large ρ, we get ψ scattering ia r→∞ = e −ipar 2ip a r (−1) +1 δ ia (2 + 1)P (cos θ) + e ip i r 2ip a r δ ia (2 + 1)P (cos θ) + 2ip a f ia (θ) .
(51) The label a represents the initial wave packet, incoming with momentum p a , while i labels the possible final states, and in general we allow for i = a. Exploiting rotational symmetry we can write down the full solution as where we have used the fact that m = 0. Then we need to find the correct asymptotic behavior that matches eq. (51). The asymptotic behaviour of the reduced radial functions that matches eq.(50) is given by The boundary conditions at the origin are such that the regular solutions are chosen u ia ∝ r +1 . Notice that the asymptotic matching (53) can be realised by requiring that at (numerical) infinity r ∞ often for the off-shell mode i = a, the condition u ia (r ∞ ) = 0 is more stable. These two conditions allows us to compute the full wavefunction for a given initial state potential. We can then extract other information, such as Sommerfeld factors and phase shifts.
Sommerfeld factors for any . The Sommerfeld factor for generic can be extracted from appropriate number of derivatives of the reduced wavefunction (recall that the Sommerfeld factor is defined as the ratio of -wave part of the wavefunction at the origin and the unperturbed wave along z): In the above expression we have made use of the fact that spherical Bessel functions expand as j (ρ) → ρ /(2 + 1)!! for small ρ, and exploited the knowledge that u ∼ r +1 at the origin. For the specific case of s and p-wave scattering, we get Phase shifts for any . The phase shift of the incoming wave packet can be computed from the amplitude of the out-going function. From eq. (50) it follows that, which in turn allows us to compute f (θ) as The cross section for elastic scattering then given by, C Modelling the nuclear potential The computation of bound state formation with wavefunctions requires the choice of a potential. In regime of shallow bound states the result is only weakly sensitive to details of the potential so one can choose the most convenient. As from textbooks on quantum mechanics [37], a convenient parametrization is the spherical well potential, with tunable parameters. Since the strong interactions are isospin symmetric, in each channel of spin and weak isospin we have a corresponding potential V (r) = −V N θ(r 0 − r), with V N > 0 and where r 0 is the range of the interaction, that we assume to be related to the dark pion mass by r 0 ∼ M −1 π . The only difference among different channels will be the depth of the well. This is the potential that is included in the numerical simulation. By knowing the spectrum of dark bound states in a given channel of spin and isospin, is then simple to include this effects in the initial state wavefunctions, by simply choosing just V N to reproduce a given binding energy. For a fixed range r 0 is always possible to tune the depth of the well V N in such a way to have only one bound state with zero angular momentum and arbitrarily small binding energy. The normalized reduced wavefunction with binding energy E B = γ 2 /M is explicitly given by while the binding energy is given implicitly by the solution of When the potential supports a single shallow bound state the equation above can be solved as, where x 0 = 1/(M r 0 ). Some benchmarks values are, In our numerical simulation we fix x 0 = 1 corresponding to large pion masses M π ∼ M and E B /M in a range between 0.1 and 0.01.
s-wave process. For s-wave the normalized positive energy wavefunction is The phase shift is determined by regularity of the wavefunction at r = r 0 , In the limit of small binding energy using (61) one finds This result can be directly derived using the effective range expansion. The amplitude for elastic scattering has the general form In the low velocity regime the scattering phase admit the expansion p 2 +1 cot δ = −a −2 −1 l + O(p 2 ) where a is the scattering length. Thus the amplitude has a pole for i √ M E B = 1/a that can be trusted for E B M and eq. (65) follows. From the wavefunction above we can extract the Sommerfeld factor using eq. (55). The result is where we have used p cot δ 0 ≈ −1/a 0 effectively resumming a 0 p to all orders. For vanishing binding energy the SE diverges in the low velocity regime giving rise to the peaks in the cross section. For magnetic transition the overlap integral between initial and final state is given by, where in the last step we expanded in the low velocity limit. Substituting in eq. (17) this gives the magnetic cross section in eq. (20). Note that to leading order the result is independent of r 0 so that the computation can be simply performed in the limit r 0 = 0 [27].
p-wave process. Let us repeat the exercise for = 1. The condition for a zero energy bound state is now The explicit solution is u p (r) = N sin(Kr) Kr −cos(Kr) θ(r 0 −r)+ √ 12π p cos δ 1 sin(pr) pr −cos(pr) −sin δ 1 cos(pr) pr +sin(pr) θ(r−r 0 ) where N and δ 1 are determined by matching the wavefunctions and the first derivative at the boundary. The SE amplitude is then For shallow bound states we can get an explicit formula as follows. To leading order tan δ 1 = −(pa 1 ) 3 . By looking at the elastic amplitude a ≈ 1/ √ M E B . Plugging in the formula above we find Note the very different energy dependence compared to s-wave processes. We can compute the cross section for the production of an s-wave shallow bound state from a pwave initial state using the explicit wavefunctions (70) and (63). The relevant matrix element is given by that does not depend on the initial scattering length to leading order.

D Modified Variable Phase Method
(Based on procedure in [38] as modified by [39], with modified boundary condition from [24].) We want to solve the radial Schrodinger equation for the reduced wavefunction u , [u (x)] in + 1 + ( + 1) where x is a dimensionless radial variable defined as x = M βr, and i, n = 1, · · · , N corresponding to the N different components of the wavefunction. For the EW triplet case, N = 2 for the charged-charged and neutral-neutral components of the (spin-0, charge-0) wavefunction. We split up the potential into an asymptotic part and a short-distance part as We will use the regular and irregular solutions (f i (x) and g i (x), respectively) of the 'free' Schrödinger equation: The free solutions are normalized such that the Wronskian, i.e.
for dimensionless wavenumberk i = (1 − V ∞ ii /(M β 2 )) 1/2 . The variable phase ansatz states that we can write the regular solution to the full Schödinger equation as linear combinations of the 'free' solutions for unknown functions α in (x), β in (x) for each , with boundary conditions α in (0) = δ in and β in (0) = 0. We have doubled the degrees of freedom, so we also impose the following constraint which reduces the N second-order equations, Eq. (74) to a system of 2N first-order coupled ODEs: and β(x) is linearly dependent on α(x): Using (80) in (79) we obtain: where the argument x has been suppressed everywhere for brevity. Introducing matrix N ij minimizes numerical convergence issues (see [38] for more details): Differentiating the above, and substituting for O from (81), we can derive that N ij satisfies the following differential equation: Following [38] we can write [u ] in = N j=1 N ijα jn g n whereα in = gn g i α in . Differentiating this, we obtain a differential equation forα in that is logarithmically stable:α So to solve for the physical wavefunctions, we need to solve first-order differential equations (82) and (83), with boundary conditions (84). Unfortunately imposingα boundary conditions at small x leads to numerical instabilities, and so we impose instead lim x→∞α in (x) = δ in .