Cosmological production of dark nuclei

We study the formation of Dark Matter nuclei in scenarios where DM particles are baryons of a new confining gauge force. The dark nucleosynthesis is analogous to the formation of light elements in the SM and requires as a first step the formation of dark deuterium. We compute this process from first principles, using the formalism of pion-less effective theory for nucleon-nucleon interactions. This controlled effective field theory expansion allows us to systematically compute the cross sections for generic SM representations under the assumption of shallow bound states. In the context of vector-like confinement models we find that, for nucleon masses in the TeV range, baryonic DM made of electro-weak constituents can form a significant fraction of dark deuterium and a much smaller fraction of dark tritium. Formation of dark nuclei can also lead to monochromatic photon lines in indirect detection. Models with singlets do not undergo dark nucleosynthesis unless a dark photon is added to the theory.


Introduction
The stability of Dark Matter (DM) on cosmological time scales strongly suggests the existence of new accidental symmetries in Nature. In a minimalistic approach where DM is a representation of the Standard Model (SM) gauge group, the possibilities that DM is accidentally stable are very limited and constrained [1]. A different class of models where DM is accidentally stable with a wide variety of SM quantum numbers is the one where DM is a baryon-like object of some strongly interacting dark sector [2][3][4], in the framework of vector-like confinement [5]. In these scenarios, the presence of elementary 'dark' fermions with an accidental U(1) symmetry guarantees the stability of the lightest dark baryons, as the stability of protons and nuclei is related to baryon number conservation in the SM.
In such scenarios, dark nuclear forces are expected to give rise to dark nuclei. For example in [6,7] it was shown that bound states with baryon number 2 exist and the absence of a Coulomb barrier implies that states with large baryon number very likely JHEP04(2019)108 exist in the spectrum. The possibility that DM is confined into more complex structures is of obvious theoretical and experimental interest.
In this work we wish to quantitatively study the nucleosynthesis of the dark sector. In the SM, the formation of light elements during Big Bang nucleosynthesis (BBN) depends upon a few odd circumstances: for example on the fact that the deuterium binding energy is small and comparable to the proton-neutron mass difference. We find that even for DM made of dark baryons, the success of dark nucleosynthesis depends on the formation of dark deuterium but, contrary to the SM, the different densities and binding energy of DM require a precise knowledge of the dark deuterium production cross-section.
First, we establish the general features of dark nucleons and nuclei in models with a new confining gauge interaction and fermions in vector-like representations. Under broad assumptions the spectrum is dominated by the nuclear binding energies E B much smaller than the nucleon mass M N and electro-weak effects can be included perturbatively. Importantly the synthesis of nuclei requires release of energy which automatically allowed in models with electro-weak constituents where it is always possible to radiate a photon, The key input to determine the abundance of dark nuclei is the dark deuterium crosssection. At first sight its calculation involves strongly coupled nuclear reactions that seem difficult to control. This is not the case, however, in light of the smallness of E B /M N and a precise computation is possible for shallow bound states. We will get inspiration from effective field theory of nucleon interactions [8], applying the same techniques to the case of dark nucleons with arbitrary quantum numbers (r, S). The cross-section for bound state formation through electric (dipole) and magnetic interactions are found, where the first factor accounts for possible Sommerfeld enhancement (SE) and K E and K M are group theory factors. With these result in hand we can easily compute the abundance of dark deuterium and heavier dark nuclei in a given model by solving the relevant Boltzmann equations. For models with electro-weak charges we find that a fraction of DM is bound in dark deuterium and heavier dark elements are not significantly produced. Formation of dark nuclei through photon emission can be tested in indirect detection experiments such as FERMI.
The present work clarifies how in general the first steps of dark nuclei formation can occur, based on first principle calculations. In this respect it provides a quantitative input for works that focus on the asymptotic fusion of large dark nuclei [9,10], and also for models where the dark baryon coupled to dark photons as in [11], that we reconsider. On a JHEP04(2019)108 more technical side, this paper extends the analysis of cosmological bound state formation of perturbative bound states [12,13] to strongly coupled nuclei.
The paper is organised as follows. After an introduction to the properties of nuclei in vector-like confinement scenarios in section 2, we derive in section 3 the Boltzmann equations for dark deuterium in a few models. We compute the dark deuterium formation cross section with non-relativistic effective field theory techniques in section 4 and the appendix. Section 5 discusses the case of dark nuclei made of SU(2) L triplets, while in section 6 we consider dark nuclei charged under a dark photon. We conclude and outline future directions in section 7.

Dark force
We frame our discussion of DM nuclei within the scenarios of vector-like confinement [5]. The SM is extended with a new non abelian gauge force and fermions charged under the SM and dark group, described by the renormalisable lagrangian, Such framework has special interest for DM as it automatically generates accidentally stable DM candidates [2]. We will focus on SU(N) DC gauge theories with fermions in the fundamental of dark colour 1 and masses m Ψ smaller than the confinement scale Λ. Upon confinement the spectrum consists of dark pions and dark baryons with charges under the SM determined by their constituents. The theory features an accidental U(1) symmetry, the dark baryon number, under which the fermions Ψ transform with the same phase. This symmetry guarantees the stability of the lightest baryon, which for appropriate choices of SM representation can be a viable DM candidate. The very same dark baryon number also guarantees the exact stability of the lightest states in each charge sector. This leads to the stability or metastability of nuclei or more in general state of matter that carry baryon number. The quantum numbers under the SM are uniquely fixed as gauge interactions naturally select the smallest representation as the most bound.
As an example we will consider the simplest models with SU(3) DC gauge group: • Ψ is a triplet under a SU(2) L . The lightest baryon V ≡ ΨΨΨ is also a triplet and has spin 1/2. The lightest states are dark pions transforming as a triplet and quintuplet with mass splitting ∆M 2 Π ∼ α 2 /(4π)Λ 2 . The triplet is accidentally stable due to G-parity [15], but it can decay through dimension-5 operators.
• Ψ is a singlet of the SM. In this context some heavier unspecified charge fermions are needed to guarantee thermal contact with the SM but play no other role in the dynamics. The lightest baryon has spin 3/2 for one flavour of singlets or 1/2 for more flavours.

JHEP04(2019)108
Previous studies focussed on symmetric DM where the relic abundance is generated through thermal decoupling, which leads to a baryon mass around 100 TeV to reproduce the known DM density. Here we will focus on asymmetric DM as this maximises the formation of nuclei. This requires a suppressed symmetric component so naturally M N < 100 TeV. Direct searches at the LHC place a bound of around 1 TeV for the mass [16,17], while the singlets just need to be heavier than about 10 MeV to avoid bound on number of species.

Properties of dark nuclei
At zero temperature dark nucleons are expected to bind into larger nuclei due to residual strong interactions. Based on nuclear physics examples, as well as lattice results [18], we consider as typical binding energies, In the SM due to the electro-static repulsion, only nuclei with atomic number A 100 are long lived or cosmologically stable. Moreover the presence of bottlenecks associated to the details of the nuclear spectrum implies that only the lightest nuclei are synthesised cosmologically.
For accidental DM models with electro-weak constituents the situation is likely different. If the mass of the dark baryon is larger than the electro-weak splittings, it is possible to exploit the approximate SU(2) L symmetry to classify nuclei into SM representations at least for the nuclei with small atomic number A. Actually, in the limit where we neglect SM interactions the theory has a larger accidental global symmetry SU(N F ) corresponding to the dimension of the lightest nucleon representation. All the states can thus be classified into multiplets of SU(N F ). Electro-weak interactions break this symmetry splitting the lightest nucleon multiplet into SM multiplets by an amount, where I N is the isospin of the nucleon. The electro-weak splitting of nucleons induces a splitting between the dark nuclei. Since the shift above can be larger than the nuclear binding energies, the splitting of nuclei made of different nucleon representation is dominated by the splitting of constituents over a large range of parameters. One consequence is that as in the SM heavier baryons do not participate to the dark nucleosynthesis. Cosmologically the bound state formation begins at temperatures of order E B /20, where E B is the binding energy. Since the heavier nucleons decay into the lighter ones through strong interactions we can neglect the population of heavier nucleons as long as, This is the typical range expected for nuclear binding energies. Thus we can focus on the nuclei made of the lightest SM rep. The nuclei made of the lightest nucleon multiplet can be decomposed into SU(2) L reps split by nuclear binding JHEP04(2019)108 energies ∆E N B . The electro-weak correction to the binding energy is given by, where λ = I N (I N + 1) − I B (I B + 1)/2 is the effective coupling in the isospin channel of the nucleus. For shallow bound states we estimate the size of the nucleus with the scattering length 1/a = √ M N E B , see section 4. It follows that the nuclear binding energies dominate for E B /M N > 10 −3 . Electro-weak corrections can be included perturbatively in the relevant region of parameters, and for R nucleus < M −1 Z the SM gauge fields can be treated as massless.
For large A a multitude of SM reps exist with isospin up to AI where I is the nucleon isospin. Of these the smaller representations are attractive making the nuclei more bound while the largest representation have a Coulombian energy that scales as AI(AI + 1) that unbind nuclei of arbitrarily large charges. In light of this, the valley of stability of dark nuclei will likely extend to very large to A, at least for small electro-weak charges.
Only the lightest baryon in each baryonic number sub-sector will be stable so that at late times all baryons produced cosmologically decay to a neutral state with isospin 0 or 1/2. This process is controlled by the decay rates among and within isospin multiplets induced by electro-weak interactions, analogously to de-excitation of hydrogen atom. In particular, • Each isospin multiplet can decay to states with smaller isospin. Since we focus on the s-wave bound states the rate is dominated by emission of a photon through magnetic dipoles with ∆S = 1. As we show in appendix B the rate can be computed model independently as, where B 1,2 are the binding energies of two bound states with ∆I = ∆S = 1.
• The splitting within each baryonic electro-weak multiplet is given by [1], 2 which for an SU(2) L triplet gives ∆M N = 165 MeV. Since charged and neutral components remain in equilibrium, the abundance of charged nuclei is then approximately In what follows we will study the cosmological synthesis of dark nuclei. Cosmologically states with large A could be produced through fusion processes via aggregation of heavy elements. In [9,11] it was argued that at least for light DM the dark synthesis is very efficient and one ends up with a distribution of large nuclear states. These studies overlook the first step of formation of nuclei with baryon number 2 (dark deuterium) that cannot take place through simple fusion processes but requires some energy to emitted. As we will show the dark deuterium abundance is often suppressed leading to a small abundance of larger dark nuclei.

JHEP04(2019)108 3 Dark deuterium abundance
As in the SM, we assume a separation of scales between the nucleon masses M N and the nuclear binding energies E B . In this case the treatment of dark nucleosynthesis becomes relatively simple. At temperatures below E B , when dark nuclear reactions can form dark nuclei, the dark nucleons are non-relativistic and already decoupled from the SM plasma. The yield of dark matter number is then set by where n DM is the DM number density and s = (2π 2 /45)g T 3 the entropy of the plasma and g the relativistic degrees of freedom of the plasma at a temperature T . We will assume for simplicity that DM is asymmetric, since this maximises the production of bound states. Given that thermal abundance indicates M N ∼ 100 TeV we assume the DM mass to be below this value. This is also necessary to produce a significant fraction of nuclei.
Dark BBN takes place during a period where the dark baryon number is conserved, where we have introduced the mass fractions X A i = AY A i /Y DM of a nucleus with atomic number A i . Practically, we are only tracking how the total DM yield in eq. (3.1) gets redistributed among different nuclear species. A successful dark nucleosynthesis of states with A 2 depends on the efficiency of dark deuterium formation, similarly to the SM case. In this section, we derive the parametric dependence of the dark deuterium abundance on the relevant parameters of the theory with a focus on two scenarios • SM charged constituents, with production N + N → D + X, with X = W, Z, γ.
• SM neutral constituents, with production N + N + N → D + N .
In general bound state formation requires some energy to be released. We do not consider the possibility of emitting a pion of the strong sector. At first sight this process could be favoured by the strong coupling to nucleons but it is likely forbidden kinematically. Indeed, based on nuclear physics examples [18], we expect the following scaling This is particularly significant for models with singlets where no other light states exist.

SM charged constituents
When the constituents of DM have electro-weak charges, direct searches at LHC imply that the scale of new states must be larger than about 1 TeV [16,17]. 3 For this value of the JHEP04(2019)108 mass the numerical density in eq. (3.1) is much smaller than the yield of baryonic matter, suggesting that the formation of bound states less likely than in the SM. We consider the first step of formation of dark deuteron N + N → D + X, where X stands for an electro-weak gauge boson γ/Z/W in equilibrium with the SM bath. For simplicity we assume that no other channels contribute, so that the Boltzmann equation for D takes the formṅ where n B and n D are the numerical densities of dark baryons and dark deuterium and we assume radiation domination. During radiation domination, the Hubble parameter is The second term in the square brackets of eq. (3.4) becomes exponentially suppressed as e −E B /T , when the temperature drops below the binding energy. At this stage, dark deuterium can be formed, so that it is convenient to introduce a new time variable z ≡ B/T . In terms of the mass fractions, the Boltzmann equation can be cast in the following form, (3.6) where c = 1.32, β = 45/(16π 7/2 ) and g * is the effective number of degrees of freedom. Only at a temperature significantly lower than the binding energy the bound state can be produced without being immediately dissociated. Given the smallness of Y DM , this happens at values of z f ≈ 20, when the coefficient of the term linear in X D inside the square brackets becomes of O(1). At this time we have the most efficient stage of deuterium synthesis starting with boundary condition X D (z f ) ≈ 0. Soon after, X D approaches a constant value that is linearly sensitive to the product of the binding energy and the cross-section, given by Away from saturation, X D 1, we have the following for an s-wave cross section, Contrary to the SM the production of dark deuterium is far from saturation for heavy dark matter. For this reason the actual abundance depends on the precise value of the cross-section that we will compute in section 4.

Neutral constituents
When the fundamental fermions are SM singlets, the hadrons and mesons of the dark sector can be much lighter than the electroweak or even QCD scale, without occurring in JHEP04(2019)108 constraints from LHC direct searches. As a consequence DM numerical densities comparable or even larger than visible matter (see eq. (3.1)) become possible. Naively this is the most favourable situation for dark nucleosynthesis and even very large dark nuclei could be formed as advocated in [9].
Nevertheless in this section we would like to argue that, in presence of SM singlets, dark nucleosynthesis is very unlikely to take place unless extra light degrees of freedom such as a dark photon [11,20] or a scalar [21][22][23][24] are included (see also realisations in scenarios of mirror world [25]). In absence of light fields external to the strong sector, such as a dark photon or light dark pions (see the discussion in the previous section), the first step of dark nucleosynthesis cannot occur as a 2 → 2 process. The dark deuterium production must necessarily proceed through 3 → 2 processes involving only baryon states, such as N + N + N ↔ D + N reactions. Even at temperatures below E B , when the production reaction is the only one that can occur, the fusion is suppressed as compared to the previous case by an additional power of Y DM . The Boltzmann equation in this case takes the forṁ where we have introduced a generalised 3 → 2 cross-section with mass dimension −5. The thermally averaged 3 → 2 cross-section is defined as (3.10) For the analog of s-wave processes at low energy σv 2 goes to a constant. Following [26] we estimate Introducing the deuterium baryonic fraction the Boltzmann equation (3.9) can be written as where a = 1.16 is a numerical coefficient arising from the evaluation of H and s. We can now compare the third line of the above equation with eq. (3.6) that sets the abundance of dark deuterium from 2 → 2 processes. Compared to section 3.1 for z z f the source term decouples as fast as z −4 , is suppressed by higher powers of the binding energy and especially by one more power of Y DM .
We then obtain the estimate X D ∼ 10 −13 g 10 which is utterly negligible in the relevant range of parameters.

JHEP04(2019)108 4 Production cross section of dark deuterium
In this section we explain how the cross-section for formation of nuclei can be computed from first principles exploiting universal properties of short range nuclear interactions. The main point is that, for shallow bound states such as nuclei where E B M N , it possible to write a general effective theory of nucleons. This effective field theory (EFT) reproduces the effective range expansion of quantum mechanics and allows us to compute the properties of nuclei such as their production cross-section, see [27] for a review. Here, we quickly outline the formalism and apply to the formation cross-section of dark deuterium.

Dark nucleon effective field theory
The production of nuclei with A = 2 is a process entirely analogous to the deuteron formation in the SM, pn → dγ. This can be calculated in quantum mechanics with appropriate potentials, but, as noticed long ago [28,29], it does not depend much on the details of the potentials used, as long as they are short range. As emphasised in [30,31], the generality of this phenomenon is immediately captured by the π-less effective theory of non-relativistic nucleons [8] that we briefly review. We refer the reader to the appendix and refs. for more details. 4 Since the energy scale relevant for nuclei formation is much below the pion mass it is useful to describe nucleons with a non relativistic lagrangian where the dark pions are integrated out, the π-less EFT [8]. Such theory is extremely simple because it only contains contact interactions among nucleons and couplings to SM gauge fields. Generalising the results of the refs. above, the nucleons, in a generic isospin representation r, are described by the effective lagrangian where the covariant derivative D µ = ∂ µ − ig 2 A a µ J a contains the minimal coupling to SM gauge fields and we have included 4-nucleons interaction and a magnetic dipole 5 interaction where J a is generator of SU(2) L in the nucleon rep. The coefficient κ is expected to be of order unity (in the SM the isovector nuclear magnetic moment is κ V = 2.35).
At sufficiently low energies the leading interactions in L 4 include only operators without derivatives, that can be decomposed into spin and isospin channels as where the matrices CG r and P S act on the isospin and spin space respectively. 6 Remarkably, the lagrangian above also describes the non-perturbative bound state allowing to compute for example the production cross-section. A quick way to derive the 4 In [32][33][34] this formalism was applied to Wino scattering and annihilation. 5 In the case where the strong sector violates CP through a θ angle an electric dipole is also present producing similar effects for deuterium formation. We will neglect this term. 6 For SU(2) the matrices CG can be identified with the Clebsch-Gordan coefficients, while the explicit JHEP04(2019)108 main result is the following. The non-relativistic amplitude for the elastic scattering of two nucleon has the general form in each isospin/spin channel r, S, where δ r,S is the phase shift and p = √ EM N is the nucleon momentum in the center of mass frame. For s-wave scattering in the low velocity regime one can show that p cot δ r,S = −1/a r,S + O(p 2 ), where a r,S is the scattering length. This is known as the effective range expansion. When a r,S is large and positive it follows that the amplitude has a pole at negative energy. From this we can recover the general relation between the scattering length of the binding energy of shallow bound states, where E B is the binding energy of lightest s-wave bound state. The coefficients of the 4-Fermi interactions in eq. (4.2) must be fixed to reproduce the effective range expansion of nucleon-nucleon elastic scattering. As shown in the appendix, to leading order in the derivative expansion but to all order in the scattering length a r,S , on finds Once this matching has been performed, the above lagrangian can be used to compute other processes, such as the production of deuterium (see [30] and refs. therein for the SM case). Indeed, the amplitude above determines also the coupling of two nucleons and deuterium as where γ = 1/a. This effective coupling can then be used to compute the interaction between nucleons and the deuterium. Amazingly, we just need to study the elastic nucleon-nucleon (N N ) scattering to infer all the quantities needed to perform the leading order calculation of the deuteron formation rate. Since this process occurs cosmologically at energy much below the mass of the dark pions (see also the discussion in section 2.1), it is very reasonable to use an effective field theory of non-relativistic dark nucleons (with SM quantum numbers) without the dark pions. expression of the spin projectors onto the spin 0 and 1 states are The labels r and S identify the SU(2)L and spin representations, while M and i are the indices of such representations.

Magnetic and electric transitions
The effective lagrangian in eq. (4.1) together with eq. (4.7) allows us to compute the short distance cross section for the process N + N → D + V a in terms of the binding energies and scattering lengths. The nucleus can be formed through emission of a SM gauge boson either through the electric coupling (minimal interaction) or through the magnetic dipole interaction in eq. (4.1). For κ ∼ 1, as expected for strongly coupled baryonic states, the two processes can have similar size. Importantly, different selection rules apply to electric and magnetic transition so that ∆L = 1 for the first and ∆S = 1 for the second. This implies a different velocity scaling of the cross-sections.
The amplitude for the formation of bound state can be simply computed with Feynman diagrams of the non-relativistic effective theory (4.1) using eq. (4.1) for the overlap of final state with the deuteron. 7 The only subtlety arise for large scattering length of the initial state where an extra long distance contribution must be taken into account enhancing the tree level magnetic crosssection. This effect is discussed in detail in appendix A. The amplitudes for bound formation can be conveniently decomposed in the basis of total spin and isospin of initial and final states (N N, D) using the projectors of eq. (4.2). In the limit v rel 1 one finds 8 (4.10) M, M are the indices of the isospin representations, while i, i are the indices of the total spin representations of initial and final states, and so these expressions should be read taking into account the selection rule of the magnetic and electric transition. The group theory factor is [13] C aM M (4.11) The formulae above are general, and can be even applied to other gauge groups and different representations for the constituents of the bound state with minor modifications. For any more details we refer the reader to the appendix A. 7 For Majorana bound states a factor √ 2 must be included in the amplitude to account for the normalisation of the wave-function of bound states made of identical particles. 8 We neglect here the effect of non-abelian interactions [13]. This is justified if the nuclei are dominated by the strong interactions so that their size is smaller than the Coulombian Bohr radius a −1 0 = λα2MN /2.

JHEP04(2019)108
With our normalisations, the cross-section for the bound state formation can be computed with This gives the following magnetic and electric cross-sections: Magnetic cross-section. The averaged cross-section for the production of an s-wave bound state with energy E B and isospin quantum number (I , M ) through the emission of a (massive) vector boson V a from an initial state with (I, M ) at low velocity is given by, Here g N = 2(4)d R is the number of degrees of freedom of the nucleon initial state for Majorana (Dirac) particles. If the initial state supports a weakly bound state with a i = E B i M otherwise a i is negative and can be large. For example in the SM the second term is a i γ f ≈ 5 so it dominates [30,31].
At low velocities σv rel is constant as this correspond to a an s-wave capture with ∆L = 0 and ∆S = 1. The presence of a coherent contribution from the initial state scattering length can be understood by noticing that being a s-wave process, the initial state can be in principle have an unnaturally large coefficient in eq. (4.1) that need to be kept into account to all orders.
Electric cross-section (as in dipole approximation). The averaged cross-section for the formation of an s-wave shallow bound state through dipole emission of a photon is found to be, The velocity suppression follows from the fact that in dipole approximation ∆L = 1 so that an s-wave bound state is produce from a p-wave. Note that the formula above differs by a numerical factor from the cross-section for the formation of Coulombian bound state [13]. This is because the energy levels of 1/r potentials cannot be treated as shallow bound states. From the formulae in each isospin channel for electric and magnetic transitions we can recover the component cross-section, relevant for example for indirect detection in 5.3 using the appropriate Clebsch-Gordan coefficients. Let us note that the formulae above are actually general and also apply to bound state made of different representations and different global symmetry group.

SE(s)
The discussion so far has neglected the long distance Sommerfeld effect due to electro-weak interactions on the initial state. In a given SU(2) L channel, the cross section is multiplied JHEP04(2019)108 by the corresponding factor depending on whether we deal with s or p-wave initial states, see for example [35]. For massless mediators one finds, Here α eff is the effective strength of the electro-weak forces in a given channel. Importantly, taking into account this effect, both the magnetic and electric transition rates have the same scaling with velocity, σv rel ∝ 1/v. Note that such enhancement of electric transitions is not present for deuteron in the SM so that the magnetic transition dominates at very low velocities. This can be different for dark nuclei. The approximate scalings above are accurate for v rel α. The effectiveness of the SE, then crucially depends on the typical velocity during dark nucleosynthesis. Dark deuterium forms at T ∼ E B /z f where z f ∼ 20, which implies v rel ∼ E B /M N /5 at that time. Numerically the enhancement is more pronounced when the interaction are SU(2) L symmetric, as in this case the relevant coupling is α 2 , instead of α em .
Let us now discuss domain of validity of the massless approximations. From the point of view of the initial state particles, the mediator masses can be neglected as long as the de-Broglie wave-length is smaller than the range of the electro-weak interactions M −1 W . Therefore Sommerfeld effects are maximal if the above two conditions are met For v rel M W /M N the vector boson masses must be taken into account. These has 2 effects, the first to freeze the enhancement to a constant value corresponding to the critical velocity, This approximation gives results comparable to the analytic formulas derived using the Hulthen potential [13]. In addition the mediator mass produces peaks in the cross-section. As well known the resonant behaviour originates from bound states of zero energy supported by the potential at the critical mass. Around the peak the SE has the model independent form [36], where V 0 is the typical energy and E B the energy of the bound state close to threshold in the initial state. We also note that in the limit of large scattering length of the initial state the second term in (4.10) can be interpreted as SE due to the strong interactions. Indeed this is large when the initial state supports a bound state with energy E B M N . Using 1/a = √ E B M N and extending the computation of bound state formation to finite velocity of dark nucleons one can check that the formula for bound state production through magnetic JHEP04(2019)108 interaction (s-wave) reduces to eq. (4.19) where E B is the energy of the bound state in the initial state while V 0 is the binding energy of the produced bound state. The subleading terms in eq. (4.13) is associated to the failure of factorisation between short and long distance effects.

A case study: SU(2) L triplet
In this section we compute explicitly the production of deuterium in the V model where the constituents are triplet of SU(2) L . In this scenario, in the limit of vanishing SM couplings, the models enjoys an SU(3) F flavour symmetry and the lightest baryon multiplet is an octet. This decomposes under SU(2) L as where with abuse of notation we named the triplet nucleon V . SM gauge interactions split the two multiplets as in eq. (2.3). The triplet is expected to be the lightest state in light of the smaller weak charge. Since at the temperatures relevant for bound state formation the abundance of the quintuplet is exponentially suppressed we can focus on the nucleon V .
In absence of electro-weak interactions dark nuclei with dark baryon number 2 belong to the product of two octets 8 × 8 = 1 + 8 S + 8 A + 10 + 10 + 27 S , where S(A) refers to the symmetry of the isospin wave-function. Lattice studies indicate that all these representation actually form bound states [18] and also provide the binding energies in different channels. As expected in the flavor symmetric limit the singlet (corresponding to the H-baryon in QCD) is the most bound. Following the discussion of section 2.1, we can work in a limit where we can classify the bound states according to SU(2) L representations, while neglecting the baryon in the 5 whose abundance is Boltzmann suppressed. Therefore the dark deuterium of this model belongs just to the product Anti-symmetry of the full wave-function implies that singlet and quintuplet of SU(2) L (D 1 and D 5 ) are spin-0 while isospin triplet are spin-1 (D 3 ) (for s-wave bound states). The triplet and quintuplet deuterons are branches of 8 A and 8 S respectively so they are heavier than the singlet, with splitting dominated by the strong interactions. The classification can be generalised to larger dark nuclei. The dark tritium made of 3 triplets has quantum numbers, where the symmetry of the wave-function determines the spin. We can estimate the electroweak binding energy by adding a nucleon to the deuterium and taking into account the reduced mass. We summarise the bound states up to dark Tritium in table 1.

JHEP04(2019)108
Name r S λ Constituents Table 1. On the left quantum numbers of dark nuclear bound states (deuterium and tritium) for nucleons SU(2) L triplets. On the right group theory factors for the transition from an initial state with isospin r to a bound state with isospin r .

Relic abundances of deuterium made of SU(2) L triplets
To compute the abundance of deuterium we assume that the initial state is SU(2) L symmetric, i.e. we neglect the mass difference between V ± and V 0 . Using (4.13) and (4.14) the cross-section for deuterium formation through emission of W, Z, γ, averaged over initial states, is approximately given by, The first term in bracket corresponds to magnetic ∆S = 1 transitions and the second the electric one ∆L = 1. The triplet can be produced either from a singlet or quintuplet channels. The scattering length are given by 1/a i ≈ M N E B i . If a state is unbound the formulas above still apply with a negative scattering length. This is the case of deuterium in the SM where nn is weakly unbound producing a large scattering length.
In figure 1 we show the numerical values of the electric and magnetic cross-sections 1 ↔ 3 normalised to σ 0 = πα 2 /M 2 for various choices of the binding energies. Due to the SE described in section 4.3 the p-wave electric cross-section can be larger than the magnetic one.
To compute the abundance of deuterium in principle one should write a different Boltzmann equation for each bound state. We can however simplify the problem by noting that transition between different bound states are fast so that they are in equilibrium among them (see section 2 and appendix B). This implies that n D i /n D j = n eq D i /n eq D j . The abundance of deuterium is then determined by eq. (3.6) with the effective cross-section and degrees of freedom, In figure 2 we present the total mass fraction of deuterium X D as a function of the M N for different choices of binding energies of the singlet (D 1 ) and triplet (D 3 ) deuterium. We assume for simplicity that D 5 does not play a role even though being the least bound isotope with baryon number 2 it could enhance the production of D 3 through a large scattering length. Solid lines correspond to the abundance including electro-weak SE effects while dashed lines are obtained with the short distance cross-sections. The SE enhancement has significant impact especially because it eliminates the velocity suppression of electric transitions. For E B < M W only the photon can be emitted with smaller cross-section in light of the electric coupling and multiplicity. The transition between SU(2) L symmetric emission and photon emission originates the features in the plot.
Differently from the SM most of the baryons can form deuterium only for large binding energies. An order 1 fraction of deuterium can only be obtained for TeV masses, around the experimental collider bound. The large mass scale associated with DM with electroweak charges is the principal obstruction to convert into deuterium and then heavier nuclei an O(1) fraction of DM. Notice that there is no reduction in the mass fraction of heavy nuclei caused by the reduction of V ± from the decay V ± → π ± V 0 , which happens only T < 100 MeV. On the contrary in the SM the main limitation to form deuterium is the neutron decay.

Production of nuclei with A ≥ 3
The formation of dark tritium T is determined schematically by the following reactions where we allow for weak and strong T formation processes, the latter without the emission of SM radiation. As for the deuterium, dark tritium is produced via exothermic reactions when E T − E D > 0 (weak process) and/or E T − 2E D > 0 (strong process). If they hold, when the temperature drops below E T − E D and E T − 2E D the dissociation of dark tritium is exponentially suppressed. The most favourable situation arises, as in the SM, for E T − 2E D > E D , such that at the time of dark deuterium formation, the dissociation of dark tritium is already ineffective. These assumptions greatly simplify the discussion and they allow us to estimate the upper bound of dark tritium abundance. The Boltzmann equations including deuterium and tritium are given by Where we have introduced the following notation In deriving the above Boltzmann equations we have written effective production rates including the effect of nearby bound states with the same baryon number. In particular,

JHEP04(2019)108
σ D v eff is defined in eq. (5.6) and the reactions for Tritium are defined similarly, although we would like to distinguish the weak σ T v eff from the strong σ T v strong eff process. This inclusive approach allows us in principle to take into account all possible bound states of table 1.
When the dissociation rate of dark deuterium, D +W → V +V , becomes exponentially suppressed, we see that the terms proportional to b 1 and b 2 tend to transfer a fraction of X D to X T with an overall decoupling as fast as 1/z 2 (neglecting possible enhancements from Sommerfeld effect). As expected the strong fusion reactions (b 2 -term) have smaller rates per deuteron since they are proportional to X 2 D , partially reducing the increase in the hard rate, b 2 /b 1 ≈ 1/α.
Knowing the reaction rates one can simply solve numerically the above set of equations, however it is interesting to analyse it in the limit of small nuclei mass fractions (i.e. X D 1). For the most realistic cases we expect z 0 ≈ O(1), therefore the non-trivial evolution happens when the overall rate is very small z 0 /z 1. By contrast, in the BBN we have z 0 ≈ 10 4 , which would allow for a fast fusion of heavier nuclei if bottlenecks for the binding energy were absent [37]. The formation of dark Tritium is then further suppressed and we have These expression are accurate as long as one can neglect the loss terms in the Boltzmann equation for X D and at leading order in b 1 and b 2 . This is reliable as long as

Indirect detection
The possibility of forming bound states is relevant for indirect detection as it would lead to the emission of monochromatic photons of energy equal to the binding energy of the nucleus. This possibility is of great interest also because asymmetric scenarios typically do not produce indirect detection signals since the DM cannot annihilate into SM particles. It is also independent on whether dark nuclei are synthesised cosmologically. We leave a detailed study to [38] and here outline the main results. In the V model at late times DM is made of the neutral component V 0 of the nucleons plus a model dependent population of deuterium, also in the neutral component. We will neglect the deuterium population for the purpose of this section. If the nuclear binding energy of deuterium is larger than M W one can form the triplet deuterium via the tree-level The magnetic transition (4.13) including long distance nuclear effects (see below) gives (5.13) where we neglected the electric transition which is velocity suppressed in absence of SE. For E B 3 < M W the W is off-shell leading to a suppressed cross-section.
Another contribution arises from the SE due to SU(2) L gauge interactions and dark nuclear forces. Thanks to this effect, two neutral dark nucleons can form dark deuterium (in JHEP04(2019)108 neutral component) through the emission of a photon. Physically this is possible because |V 0 V 0 S is not a mass eigenstate and can oscillate into |V + V − S . In the symmetric limit the mixing can be extracted from the Clebsch-Gordan coefficients, (5.14) The following processes are then possible such that we have formation of neutral D 3 from an initial state |V 0 V 0 0 in the singlet or quintuplet of SU(2) L , either from p-wave spin-1 through an electric transition or from s-wave spin-0 via a magnetic transition.
Indirect signals from Sommerfeld of SU(2) L gauge interactions. At the low velocities relevant for indirect detection the electro-weak symmetry breaking effects can be important and a numerical solution is required, see [39,40]. For s-wave the SE due to electro-weak interactions is identical to the one of Wino DM χ studied widely studied in the literature, see [41]. In particular we are interested in the long distance effects that allow the neutral Winos state χ 0 χ 0 to oscillate into χ + χ − . We define as SE 00→+− the corresponding Sommerfeld factor of the Wino that can be derived from the following potential where A = α em /r+α 2 c 2 W e −M Z r /r, B = α 2 e −M W r /r and ∆M is the mass splitting produced by electroweak symmetry breaking, equal to ∆M = 165 MeV.
Noticing that in our case we have the same SE of Wino DM, so that we can exploit the calculation performed for the Wino and apply it to the dark deuterium. Therefore, the indirect detection signal can be simply calculated as The short-distance contribution to the cross-section of the charged components can be computed using (4.13) and (5.14),  Figure 3. Indirect detection bound from emission of photon lines of energy E D3 through magnetic coupling (κ = 1). Exclusion limits from FERMI are obtained by recasting the results in [43] to account for the different ratio of photon energy and DM mass. SE for the triplet has been extracted from [41].
The annihilation rate of Winos into photon pairs, σ χ 0 χ 0 →γγ+ 1 2 γZ , has the property that both the γγ and γZ final states are reached from χ 0 χ 0 with the same SE 00→+− . This allows us to provide an alternative form for eq. (5.18) in terms of Winos cross-sections [39] Due to SE the effective cross-section has peaks whose location is identical to the one of the Wino. The first peak appears for M N ≈ 2.3 TeV, but the energy of the emitted photon is now equal to the binding energy of the bound states rather than DM mass. The sensitivity to these photons can be studied according to the recipes described in [42]. In figure 3 we recast the bounds from FERMI from the galactic center [43]. Formation of deuterium through magnetic emission of photons is strongly constrained for E B > M N /10.
Indirect signals from strong interactions. Beside the standard SE due to electroweak interactions also the strong interactions can enhance the cross-section when the initial channel supports a bound state at threshold. This effect is captured by the scattering length in eq. (4.13). We can compute the cross-section relevant for indirect detection in the SU(2) L symmetric limit. Since the quintuplet is expected to be less attractive or even weakly repulsive the largest effect will be associated to the quintuplet initial state. Using the Clebsch-Gordan coefficients of eq. (5.14) and neglecting further enhancement from electro-weak interactions we find, Strong interactions favour the SU(2) F singlet as the lightest state (the analog of the SM deuterium) that will then be absolutely stable. Differently from the SM also the triplet could be bound. We will assume that the temperature of the dark sector is of the order of the SM bath which can be realised introducing heavy fermions charged under the SM. Singlets allow the mass scale to be much lower than TeV. As discussed in section 3.2 despite the larger numerical density of nucleons deuterium cannot be formed because 3 → 2 processes are very suppressed. This conclusion can be avoided introducing a dark photon that carries away the energy. Leaving the model building aspects to future work we assume the fermions to have opposite unit charges, Q D = 2J 3 so that the nucleons have charges ±1. This could be modified with different charge assignments with minor changes to the nuclei formation.
Neglecting model dependent SE due to the dark photon the cross-section for formation of deuterium through emission of a dark photon is, for the doublet rep and assumed m γ E B . In figure 4 we show the abundance of deuterium obtained integrating the Boltzmann equation and assuming that no other nuclei are produced. For masses in the GeV range, suggested by the coincidence of DM and visible matter density, the production of deuterium is very efficient even for dark photon couplings as small as 10 −4 . Note that electric transitions, even though velocity suppressed, are relevant for small binding energies as the one of deuterium in the SM (E B /M N = 0.0022). Differently from the SM for degenerate masses we do not expect significant bottlenecks so that production of heavy elements can proceed unsuppressed. For 100 GeV masses the production of deuterium is again suppressed due to the smaller numerical density of baryons.

JHEP04(2019)108 7 Conclusions
Big Bang dark nucleosynthesis is a cornerstone of the cosmological history of the Universe and one might wonder if a similar process could take place in the dark sector. In this work we have studied the synthesis of dark nuclei in theories where DM is a baryon of a new gauge interaction. Such models, motivated by the idea of the accidental stability of DM, also predict the existence of stable and metastable nuclei with SM charges that could be formed during the evolution of the Universe.
The key step for dark nucleosynthesis is the formation of dark deuterium, a nucleus with dark baryon number 2. This process requires energy to be released and it is automatically possible through emission of SM gauge bosons in models where the fundamental constituents have electro-weak charges. Remarkably, the relevant cross-section can be determined from first principles in terms of the binding energies of nuclei. To this aim, we have found it useful to employ pion-less EFT for the nucleons. This construction reproduces the effective range expansion of quantum mechanics with short range potentials [8] and it allows us to compute in general the cross-sections for the production of shallow bound states expected for nuclear interactions. We note that the cross-sections differ from the ones often used in the literature, depending in a non-trivial way on binding energies and velocity. In particular electric transitions grow for small binding energies while magnetic transitions decrease. Electric transitions moreover are velocity suppressed unless enhanced by Sommerfeld enhancement.
Having determined the relevant cross-sections one can solve the Boltzmann equations for the abundance of deuterium and heavier elements. In the case of DM with electroweak charges, for example a baryon triplet of SU(2) L here studied in detail, one finds that only a fraction of DM binds into deuterium unless binding energies much larger than in nuclear physics are assumed. The main reason for this is simply the numerical density of DM: direct searches imply that the DM mass is greater than 1 TeV for electro-weak constituents. The small numerical abundance compared to the SM nucleons suppresses the production of nuclei. Therefore dark BBN ends after deuterium formation leaving a fraction of deuterium plus small traces of tritium. In minimal models with singlets, production of dark deuterium is kinematically forbidden and dark nucleosynthesis cannot start. This conclusion changes completely including a dark photon: for DM masses in the GeV range, deuterium is formed very efficiently due to the large density and heavier nuclei can thus be formed through fusion reactions allowing to populate nuclei up to large atomic numbers.
While in this work we have focused on asymmetric scenarios where DM mass is a free parameter, our results can be extended to symmetric models, producing extra annihilation channels and nuclei-anti-nuclei. For the simplest scenarios however the critical abundance is reproduced for masses around 100 TeV and no significant fraction of nuclei is produced.
The formation of DM nuclei is interesting experimentally as it can lead to novel signatures in DM indirect detection, even in asymmetric scenarios and change the prediction for direct detection experiments due to the different composition of DM. The emission of monochromatic photons with energy equal to the nucleus binding energy is a smoking gun of dark nuclei that could be searched experimentally [38].

JHEP04(2019)108
This work extends the computation of perturbative bound state formation in [13] to strongly coupled bound states. We have provided general formulae, that could be used in other contexts, for electric and magnetic interactions also taking into account important long distance effects associated to bound states close to zero energy. It would be interesting to generalise this formalism to the fusion of strongly coupled bound states as well studying perturbative bound states within this framework.

JHEP04(2019)108
Alternatively, one can solve for C r,S requiring a pole in the amplitude for a given binding energy.
Expanding the amplitude at the pole we can also determine the interaction of the bound state D r,S with the two nucleons. Namely the coupling is the square root of residue of the amplitude at the pole. Therefore quite generally we can write as effective coupling (A.4) Once the coupling of the bound state to two nucleons is determined, one can compute with ordinary Feynman diagrams the amplitude relevant for dark deuterium formation. This formalism is particularly powerful as it allows to compute systematically the effects to higher orders and allows to resum effects associated to large scattering lenghts.
In this work we are interested in all the processes of the type where a bound state is formed through emission of a gauge boson. This implies the existence of selection rules: in the explicit case of SU(2) L , all the leading order amplitudes correspond to a ∆I = 1 transition. Instead, for spin and angular momentum there are two possibilities: i) a magnetic dipole interaction of eq. (4.1) which corresponds to a ∆S = 1 transition; ii) an electric transition (in dipole approximation) with ∆L = 1 that originates from the covariant derivative in eq. (4.1).

A.1 Details on the magnetic transition
In the SM the deuteroun is an isospin 0 and spin 1 state. At low energy the formation is dominated by the magnetic transition) from an initial spin singlet isospin triplet channel [31]. The selection rule implies that it proceeds from an s-wave initial state ( 1 S 0 ) so that σv goes to a constant for slow nucleons. Moreover this process is significantly enhanced by the large scattering length of the 1 S 0 channel.
Our scenario differs from the SM for the different group theory structure. To leading order for large scattering lengths of the initial state there are two diagrams that contribute, The crossed circle represents the deuteron coupling (A.4), while the filled circle corresponds to the insertion of the resummed amplitude for the scattering in the (r, S) channel.

JHEP04(2019)108
The loop integral appearing in the amplitude A 2 is finite giving, (A.7) The ratio of the two amplitudes for small energy is where we have used the fact that k = E B r . By virtue of this relation, the leading order term of the magnetic transition can be computed easily just focussing on A 1 . Notably, A 2 can be of the same order and might dominate over the first one.
In full generality, to leading order in the momentum expansion, the magnetic amplitude from an initial state N + N r,s to a final state with D r ,s + W a λ is given

B Bound state decays
The formalism above also allows to compute decay rates. Focussing on bound states with = 0 and spin 0, 1 the decay of excited states to lower levels occurs at leading order in α via the magnetic interaction. The decay rate is given by The amplitude for the transition between two bound states with ∆I = 1 and ∆S = 1, with the emission of a gauge boson, is given by the 1-loop graph in eq. (A.6) evaluated at the difference of binding energies between the two states,

JHEP04(2019)108
Open Access. This article is distributed under the terms of the Creative Commons Attribution License (CC-BY 4.0), which permits any use, distribution and reproduction in any medium, provided the original author(s) and source are credited.