Big Bang Synthesis of Nuclear Dark Matter

We investigate the physics of dark matter models featuring composite bound states carrying a large conserved dark"nucleon"number. The properties of sufficiently large dark nuclei may obey simple scaling laws, and we find that this scaling can determine the number distribution of nuclei resulting from Big Bang Dark Nucleosynthesis. For plausible models of asymmetric dark matter, dark nuclei of large nucleon number, e.g.>10^8, may be synthesised, with the number distribution taking one of two characteristic forms. If small-nucleon-number fusions are sufficiently fast, the distribution of dark nuclei takes on a logarithmically-peaked, universal form, independent of many details of the initial conditions and small-number interactions. In the case of a substantial bottleneck to nucleosynthesis for small dark nuclei, we find the surprising result that even larger nuclei, with size>>10^8, are often finally synthesised, again with a simple number distribution. We briefly discuss the constraints arising from the novel dark sector energetics, and the extended set of (often parametrically light) dark sector states that can occur in complete models of nuclear dark matter. The physics of the coherent enhancement of direct detection signals, the nature of the accompanying dark-sector form factors, and the possible modifications to astrophysical processes are discussed in detail in a companion paper.


Introduction
Most models of dark matter (DM) assume that it is made up of point-like particles, or at least can be treated as such in its interactions. However, this is not necessarily the case. In the Standard Model (SM) most of the SM-sector mass in the universe is in the form of composite states-atoms, containing nuclei, made up of nucleons, which are themselves made up of quarks and gluons. Particularly noteworthy is the remarkable richness of SM-sector nuclear and atomic physics resulting from just a few basic conservation laws, most importantly baryon-number and charge conservation, and a few relevant interactions, dominantly the short-range strong nuclear binding force and long-range electromagnetic interaction. The compositeness and variety of the states that result, atoms and nuclei, plays a vital role in many important physical processes in the history of the universe, for example, in galaxy formation and all of stellar physics. Given the ease with which bound states can arise, it is clearly important to consider the possibility that DM may exist in the form of composite states too.
In this work, we will investigate a particular (and as we will argue, very attractive) class of composite DM models-those featuring cosmologically stable "dark nuclei" (DN) of large dark "nucleon" number (DNN). In particular we focus on the case where, similar to SM nuclei, there is a relatively short-range strong "nuclear" binding force with a hard core repulsion, and there is an approximately-conserved quantum number, DNN, analogous to baryon number. Unlike the SM where both protons and neutrons are relevant, for reasons of minimality and simplicity we will take there to be just one kind of stable dark nucleon (or more generally, the differences between different nucleon types to be unimportant). Importantly, we will in addition assume that in the dark sector the analogue of the longrange electromagnetic interaction between protons is absent, either because there is no massless "dark photon" or because the stable dark nucleon is uncharged. 1 With these simple assumptions there exist a very broad range of stable DN states of varying DNN, and with a binding energy per dark nucleon that saturates at a finite value. This is similar to the SM where there are stable nuclei at multiple baryon numbers with approximately constant binding energy per nucleon, ∼ 8 MeV, over a range of nuclear sizes. Unlike the SM, however, there is no longer a Coulomb repulsion term, so that in the dark sector the binding energy per dark nucleon truly does saturate, at least until gravitational effects become important, and never turns over, unlike the Fe maximum in the SM, so there are stable DN up to very large DNN. In terms of a simple liquid-drop like model, analogous to the semi-empirical mass formula in the SM, we take the binding energy per dark nucleon to behave asymptotically as where k 1 is the DNN, the volume energy term is α, and the surface tension term β. 2 While the small bound states of such a model may have highly variable properties, 3 analogous to the special properties of D or 4 He, once the DN become large enough, their properties, beyond just their binding energy, eq.(1.1), can obey simple scaling laws. Most importantly certain interaction cross sections scale geometrically in the large DNN limit, at least when suitably averaged over a range of DNN to eliminate "magic number" and resonant effects which are of subdominant importance to the physics we study.
The main focus of this paper is the regime in which such scaling properties can determine the final number distribution of DM; namely, when larger states are built up via the agglomeration of smaller ones at low temperatures, as per SM Big Bang Nucleosynthesis (BBN). Like the BBN case, the simplest version of this process applies when the DM is asymmetric [7,8,, i.e., with a particle-anti-particle asymmetry that determines the final DM abundance, similarly to the baryon asymmetry in the SM sector, though the symmetric case is also interesting. In the following we assume, unless specifically stated otherwise, the DM to be asymmetric so the DNN will take positive integer values, and we denote the DN of various DNN as 1-DN, 2-DN, etc.
In the main body of this paper, Sections 2 and 3, we will show that the dark-sector analogue of SM BBN-which we call "Big-Bang Dark Nucleosynthesis" (BBDN)-can produce DN with DNN up to or even exceeding 10 8 , resulting in a very wide variety of striking changes to traditional DM phenomenology.
Specifically, we find the number distributions of DN resulting from BBDN satisfy remarkably simple forms. For example, for asymmetric DM with plausible underlying parameter values, DN with DNN 10 8 , can be synthesised, with the number distribution taking one of two characteristic scaling forms: If there is no substantial bottleneck of BBDN at small DNN the distribution of DN sizes takes on a peaked, in log space, universal nonpower-like functional form, independent of many details of the initial distribution and interactions. This solution acts as an attractor solution and we study how the distribution function of DN approaches this universal scaling form. On the other hand, in the case of a substantial bottleneck of BBDN at small DNN we find the surprising result that DN of even larger DNN, 10 8 , are often synthesised, again with a simple form of the number distribution. Such behaviours are studied both via numerical solutions of the relevant equations, and from physically motivated approximate analytical solutions.
Such states can give rise to a variety of interesting possibilities, including: • For reasonable parameter values effectively very heavy (≥ 10 8 GeV) DM can be produced by BBDN in the form of large DN (and with one of two possible, essentially model independent distribution functions). Such heavy DM is in contrast to the usual unitarity limit for point-like DM in the case of thermal freeze-out production [51], 4 and is not usually expected in asymmetric DM models, which seek to link the DM and visible sector abundances.
• Coherent enhancement of interactions: Processes that do not probe the individual constituents will have amplitude going as the DNN k, so scattering cross sections, in for example direct detection experiments, can in principle be enhanced by k 2 10 16 compared to the case of a single DM nucleon. Taking account of the reduction in the number density of such states, one still finds an effective increase in direct detection interaction rates scaling as ∼ k, so effectively reducing expected collider signals by ∼ 1/k for a given direct detection rate compared to the standard point like case. The physics of the coherent enhancement of direct detection signals, with the accompanying possibility of novel form factors from the dark sector, is discussed in detail in a companion paper [53].
• Inelastic processes at both "high energy", of order the DN binding energy differences, and parametrically lower energy, from long-wavelength collective excitations.
Examples of the high energy processes are ones that move between states of different DNNs-fusions and fissions, but there is also the possibility of the extended structure of states leading to a spectrum of excitations, as exists in SM nuclei and atoms [53].
• States with large spin: for large composite states, there is the possibility of interactions aligning the spins of many of the constituent states, leading to DN with nuclear spin ∼ k. This is of potential interest in, for example, interactions with SM nuclei in direct detection experiments, and capture in astrophysical objects.
Since our focus is the physics of BBDN, many of the above possibilities are either only briefly touched upon in this work, or not treated at all, and we leave their detailed study to a series of companion papers. Finally, before turning to our analysis of DN and BBDN we emphasise that while there are many specific models which can realise this kind of scenario, we will focus on regimes in which the behaviour can be broadly model-independent in nature.

Basics of dark nucleosynthesis
While present-day DM may be composed of large bound states, this is generically not the case in the early hot universe. At large temperatures, T , the entropy term in the free energy dominates and the chemical equilibrium distribution has almost all DM in small-number states. For small T , compared to the binding energies, the energy term dominates, and chemical equilibrium favours large bound states.
As we discuss in detail in later sections, starting from the situation at high temperatures, large DN may be assembled by an aggregation process where fusions dominate dissociations and fissions. Specifically, if interactions are not so weak as to be frozen out by the time the equilibrium shifts to favour larger states, then larger states will be built up until fusion reactions freeze out due to a combination of falling velocity and a falling number density from both Hubble expansion and the formation of the large DN themselves.
Since, for the DM masses we consider, the DM is dilute, if we are in kinetic equilibrium then the transition from being kept in equilibrium by dissociations, to fusion reactions dominating, generally happens fast enough to be only a small perturbation to the subsequent fusion process (this technical point is discussed in detail in Appendix A). If thermalising interactions are not sufficiently fast to reduce the energy of de-excitation products before they hit another DN, then these may cause dissociations, leading to very different behaviour from the fusions-only approximation (c.f. SM recombination). We will not consider this regime in the current paper, assuming instead that the de-excitation products decay or thermalise on fast enough timescales. 5 5 In principle, another exception to the statement that, late on, only fusions are important is when the excited states produced by the fusion processes de-excite by the emission of nucleon constituents, as occurs in the SM [54] (many Q-ball models also de-excite by losing charge, e.g. [55]). In general these emitted small-k DN are either quickly re-absorbed by larger DN and do not significantly alter the dynamics of the aggregation process, or they act approximately as an external bath with which the larger DN scatter. In the latter case, as long as enough of the small-small fusion processes can occur without involving nucleon

Freeze-out of fusions
Given that fusions dominate at low T , we can obtain an estimate of the maximum size of DN built up by the aggregation process by investigating when fusion reactions freeze out.
First, let us suppose that the last fusions to freeze out are those between large DN, and also that we end up with a 'peaked' mass distribution in which almost all of the mass is in ∼ A-DN. In this case, the rate of fusions for a single DN, per Hubble time, is Γ/H ∼ σv n A /H, where σv is the thermally-averaged fusion cross section, and n A is the number density of A-DN. Since we have DNN conservation, n A = n 0 /A, where n 0 is the DNN density. If the DN are non-relativistic, as required to be aggregating, and in thermal equilibrium (e.g., with an external bath of light dark-sector particles as we discuss later), then the DN velocity v A v 1 A −1/2 . For large A, saturation of the dark sector nuclear forces implies the internal mass-energy density, ρ b , of DN matter is roughly constant with size, and that fusion cross sections scale approximately geometrically, σ σ 1 A 2/3 , and so In general, the temperature, T b , of the DN bath can differ from the temperature, T , which sets the Hubble expansion rate in the radiation dominated era, and which we assume to be dominantly set by the more numerous SM sector degrees of freedom, as the dark and SM sectors may be essentially decoupled from each other. With this in mind, and using 10.75 where the parameters are normalised to SM values both for comparison and because such a parameter region is naturally motivated by asymmetric dark matter (ADM) models. Thus, in this scenario, if dissociation stops being important around T = 1 MeV, then, from eq.(2.1), the largest mass that could have been built up is (2 × 10 7 ) 6/5 m 1 ∼ 5 × 10 8 GeV, corresponding to a radius of ∼ 480 fm for the fiducial parameter values. 6 Beyond this size, the number density and velocity are too low for interactions to occur. We will see that, for reasonable parameters, this bound can be attained.
Note that, if we scale all of our dimensional parameters to increase mass scales by a factor λ (i.e. T → λT , ρ b → λ 4 ρ b , etc.), then the freeze-out DN mass scales as m fo → λ −7/5 m fo . Changing the mass scales of the constituents, e.g. by changing the confinement scale in a strongly coupled theory, may be expected to have roughly this effect. Thus, decreasing the mass scale of our constituents will tend to increase the masses that could be built up, mainly since larger geometric cross sections are available. emission, the mass fraction of small DN will be sub-dominant, though they may dominate the number density. 6 Note that this does not depend on the dark nucleon mass m1. If some scaling other than m k ∝ k between mass and DNN held, e.g., as in the case of some Q-ball models [10], then eq.(2.1) would still apply, but with A replaced by the ratio of the final mass to the mass for which σ1v1n0 was evaluated. Such models, in which binding energies were not a small fraction of the total mass, have quite different properties.
Alternatively, the last fusions to freeze out may be those between 'small' + 'large' DN. Compared to large + large fusions, while the possible number density of small DN is larger, the number of separate fusion events onto a given DN needed for the same increase in size is larger by the same factor. However, since in thermal equilibrium the velocities of small DN are larger, the overall rate of size increase is enhanced by that factor. To qualitatively understand this scenario, suppose that an order one proportion of the dark nucleons are in small k-DN, and the rest in large A-DN. Then, the rate at which an A-DN increases in size by A is given by where we used σv A+k = 1 4 δσ 1 v 1 k −1/2 A 2/3 , with δ parameterising a possible suppression of small-large cross-sections from the geometric value, e.g. from quantum reflection effects, as per SM nucleon-nucleus interactions. So, the maximum attainable DNN of the DN is corresponding to m ∼ 7×10 19 GeV (and scaling as m fo → λ −5 m fo ) and a radius of 3×10 6 fm for our fiducial parameter choices. Similarly, the rate of fusions onto an A-DN for a given k-DN is Γ ∼ σv n A ∼ n A n k σv n k , so the largest DN that can absorb the entire population of small DN also have approximately this same maximum size. Thus, as long as δ is not too small, and a large enough population of small DN exists for long enough, there is the possibility of producing much larger DN via this route, than via an aggregation process in which the DN at any given time are of approximately the same size. Obtaining significant quantities of DN this large requires that a population of small DN remains around until the end of the aggregation process, i.e. that small + small fusions are slow. If this is the case, but small + large fusion cross sections are roughly geometrical, and we produce a number density of 'seed' large DN ∼ n 0 /A fo , then, as studied in Section 3.4, we will end up with most of the DM mass in DN of size ∼ A fo . If we produce a larger seed density, the maximum size will be reduced proportionally, up to a cross-over point at which this size is lower than the freeze-out limit for large + large fusions, at which point we enter that regime. If a smaller seed density is produced, then most of the DN never gets through the 'bottleneck', and remains small, with some sub-dominant population of large DN up to m fo . Section 3.4 investigates these regimes in detail.

Bottlenecks, and comparison to BBN
Bottlenecks to nucleosynthesis, in the form of 'anomalously slow' small + small fusion rates for certain channels, are important in Standard Model BBN. There, only small nuclei are synthesised, with almost all of the nucleons ending up in H and 4 He, only small amounts in the other A ≤ 7 nuclei, and entirely negligible amounts beyond this. The B D ∼ 2 MeV binding energy of D means that, assuming only A = 1, 2 states are occupied, n D /n p ∼ η(T /m p ) 3/2 e B D /T , and since the baryon to photon ratio η 6 × 10 −10 , this only becomes 1 at T 0.06 MeV. However, slightly before this temperature, the D number density becomes large enough that 2 + 1 and other processes start occurring, and in fact their rates are 10 4 times the Hubble rate, as expected from a calculation along the lines of eq.(2.2).
The real bottleneck preventing the synthesis of large nuclei in the SM is the large binding energy per nucleon of 4 He compared to subsequent small nuclei. It is not until 12 C that the binding energy per nucleon exceeds 4 He. In quasi-equilibrium among smallnumber nuclei, this large binding energy means that 4 He dominates the mass fraction along with H, their ratio being set by the abundance of neutrons. The small number densities of A = 7 nuclei produced are nowhere near sufficient to make the rate of production of energetically favoured A ≥ 12 nuclei comparable to the Hubble rate, and so we freeze out with quasi-equilibrium abundances (for a textbook discussion see e.g., [56]). Note that this kind of bottleneck, involving a wrong-sign binding energy difference rather than merely a small right-sign difference as in the D example, gets worse rather than better with decreasing temperature. 7 In summary, SM BBN provides an example of a nucleosynthesis process where the presence of large binding energy differences among small-number states creates a bottleneck. In fact, there are two bottlenecks, the first of which, the D bottleneck, we get through well before fusions have frozen out, but not the second, post-4 He.
A bottleneck in BBDN may similarly have the result of preventing significant quantities of large DN from being formed. However, there is also the possibility, as discussed in the previous section, of a suitable bottleneck leading to the build-up of even larger DN than would otherwise have been produced, with a qualitatively different distribution of DN sizes. Section 3.4 explores this scenario in detail.

Aggregation process
Recapping, the cross-over from the high-T regime to that in which dissociations are unimportant occurs in a sufficiently short time that only DN much smaller than the freeze-out size are able to be built up in appreciable quantities during this period. After this, there will be a period of effectively fusion-only aggregation before we reach one of the limits discussed in Section 2.1. This section investigates the details of this aggregation process.
Evolving with fusions only, the Boltzmann equation for the k-DN number density, n k (t), is where σv i,j is the velocity-averaged cross section for i + j → (i + j). 8 7 In the D case, since nD increases exponentially with falling temperature, the number of A ≥ 3 states produced beforehand is insufficient to bypass the bottleneck. For this reason, it is somewhat difficult for the bypass process to be dominant for right-sign binding energy differences. This equation combines together the dilution of number densities due to Hubble expansion, the change in cross sections due to the decrease in DN velocities, and the change in the relative concentrations of k-DN. We can separate these processes out by moving to different variables. Writing the velocity-averaged cross sections as σv i,j = σ 1 v 1 K i,j , where the K i,j encode the relative rates of different fusion processes, the Boltzmann equation in terms of the dimensionless yields Y k ≡ n k /s, where s(t) is the entropy density, is Here we assume that entropy in conserved throughout the aggregation process, so s ∝ a −3 . This system of aggregation equations can be solved in terms of relative concentrations, Note that Y 0 is constant in time by nucleon number conservation, and the solution is normalised such that k k y k = 1 throughout. In words, the function w(t) describes how quickly the number distribution of DN changes, whereas the set of distributions moved through is determined by the form of the 'aggregation kernel' K i,j and the corresponding solution to the aggregation system eq.(3.3).

Scaling regime
The system of aggregation equations, eq.(3.3) for all k, is simplest to solve when the K i,j do not depend on w, i.e. when we can absorb all of the time/temperature dependence into the single quantity v 1 (t). The simplest case in which this is true is when the DN are in kinetic equilibrium with each other, at some temperature T b (t). k-DN then have mean- If fusion cross-sections scale approximately geometrically for large DN, then a kernel of the form is a good approximation for large DN. 9 In this expression we have, for simplicity, replaced the relative velocity, which is non-relativistic in the fusion regime, by an averaged velocity. Figure 1 shows, in both w and t space, numerical solutions for the aggregation dynamics eq.(3.3) with this kernel, starting from y 1 (0) = 1 initial conditions. What is immediately these dominantly de-excite via nucleon emission, this will just reduce their net forward rate by some factor, and, in simple cases, leaves the qualitative behaviour otherwise intact. 9 This specific form matches cases investigated in the mathematical literature [57]. noticeable is that the number distributions at different times have almost the same shape, but shifted relative to each other corresponding to an increase in average size. This arises because our kernel is homogeneous, K bi,bj = b λ K i,j with, here, λ = 1/6. Physically, if we scale up the DNN of all of the DN by some factor, this just corresponds to an overall scaling of the rate of fusions, and doesn't change the relative rates of different-number processes. It is known [58] that, for such kernels, there is generally a 'scaling solution' such that almost any finitely-supported initial conditions eventually converge to the form where k(w) ∝ w 1/(1−λ) is the 'characteristic size' of DN at time w, f is the kernel-dependent scaling solution, and the 1/k 2 factor ensures correct normalisation.
In our case, we can check whether the numerical solutions in Figure 1 obey this scaling behaviour by choosing a plausible form for k (e.g. for a peaked number distribution, we could take one over the total number of DN, ( k y k ) −1 ), and plotting k 2 y k against k/k, as shown in Figure 2. We see that, for y 1 (0) = 1 initial conditions, the distribution converges very quickly to such a scaling solution, which can then be used to extrapolate to larger w values. From the mathematical theory of aggregation developed in other contexts [57,59], x for x 1, and to have the form f (x) ∼ x θ e −const./ √ x for x 1, with θ some constant, both of which match sensibly onto the numerically obtained form. Also, the behaviour of k(w) with w matches, for larger w, the ∝ w 6/5 form predicted from λ = 1/6. This also holds for different choices of k, as expected for a peaked number distribution.

(In-)Dependence on initial conditions
As discussed in Section 2, we do not expect to start the fusions-only aggregation process with single-nucleon-only initial conditions-instead, we will have whatever was produced during the phase when dissociations were still important. Furthermore, fusion cross sections between small DN will probably not be well approximated by the geometrical scaling appropriate to large DN, for example, SM cross sections between small nuclei display very complicated behaviour. Since, as we discuss in the next section, there is only a finite w time available for aggregation due to Hubble expansion, the question is whether these initial conditions, and small-k behaviour, converge to the scaling solution quickly enough for appropriate measures of convergence. If the initial conditions include a subdominant large-k 'tail' for which large + large fusions are frozen out (see the right panel plot of Figure 3 as an example), then as long as the aggregation of smaller DN proceeds not much slower than the scaling behaviour, only a small proportion of the small DN fuse with those in the tail. To see this we can approximate the tail by purely A-DN, and use the fact that the rate at which a given small k-DN fuses into the tail is Γ k+A ∼ δσ 1 v 1 n A A 2/3 k −1/2 . So, writing the fraction of the DN in the tail as α = An A /n 0 , While α is small, so the small-k build-up proceeds like the scaling solution, we have k ∼ w 6/5 . Since A + A fusions are frozen out, n A is roughly constant throughout, so Thus, either δk max A, in which case the tail is subsumed into the scaling distribution, or else α is always 1. Note that this only occurs because k increases sufficiently quickly with w. As we shall see in Section 3.4, if that does not happen, then none of the DN may ever reach the scaling regime, and α max ∼ n A n 0 w 3 max can become of order 1, so all of the small DN can be absorbed by the frozen-out tail.
Alternatively, suppose that the initial conditions have some component for which similar-size fusions are not frozen out, e.g. the left panel of Figure 3. If these have the same average size as a monomer-only scaling solution after δw, then the eventual distribution will be close to the monomer-only solution after w max +δw. Since small + large fusions are faster than large + large, the 'memory' of the initial shape is erased fairly efficiently. 10 Figure 3 illustrates these behaviours numerically, with the left panel showing the convergence of non-frozen-out initial conditions to a slightly-further-along scaling distribution, and the right panel showing that a sub-dominant frozen-out tail has little effect on the scaling solution obtained. 10 Similarly, if the small-k cross sections are larger than those extrapolated down from large k, the effect is to start the large-k process slightly earlier, w → w + δw, while if the cross sections are smaller than the extrapolation then the solution interpolates between the scaling and bottlenecked regimes-see Section 3.4.

Real-time behaviour
The scaling distribution gives us the form of the y k (w) solution-to obtain the real-time distribution y k (t), we need to know the behaviour of w(t) from solving eq.(3.4). Generally, we are most interested in the t → ∞ limit, so want to find w(t → ∞). To obtain this we in turn need to know the behaviour of v 1 (t), which is simple in the case that the DN are in kinetic equilibrium throughout the aggregation process.
The simplest way for the DN to be in kinetic equilibrium is for them to be in thermal contact with a bath of lighter particles during the aggregation process. If this bath is relativistic, so T b ∝ 1/a (here ignoring possible changes in the number of species for simplicity), and taking v 1 = T b /m 1 , we obtain, assuming radiation domination, where T 0 is the temperature when w = 0, i.e. at the start of the aggregation process. Solving this, where the limit is T 0. The right panel of Figure 1 shows the solutions obtained at different t values, assuming this relationship between w and T , illustrating the convergence towards the w = w max solution.
Alternatively, if the hidden sector bath is non-relativistic, so T b ∝ a −2 , we obtain w(T ) 1 2 Since k ∼ w 6/5 , we have k(t → ∞) ∼ n 0 σ 1 v 1 H T 0 6/5 in both cases. This agrees with the approximate freeze-out calculation of eq.(2.1), which had Γ/H ∼ n 0 σ 1 v 1 H A −5/6 . Generally, since from eq.(3.4) dw d log a = n 0 σ 1 v 1 H , and during the radiation era we have H ∝ a −2 , n 0 ∝ a −3 , and writing v 1 ∝ a −γ , we have dw d log a ∼ a −(1+γ) , so, ∆w ∼ ∆ a −(1+γ) as we obtained in eqs.(3.10) and (3.11). Thus the bulk of the aggregation process takes of order a Hubble time to complete. This is illustrated in the right panel of Figure 1, which shows solutions at half-e-folding-time intervals. Such behaviour just comes from the freezeout properties of the interactions, so in the bottlenecked regime we take at most this long as well-in fact, as we shall see in the next section, that process may take much less than a Hubble time.

Bottlenecked regime
As illustrated in Section 3.1, if the fusion rates, as parameterised by w max and the K i,j , are not too small compared to the support of our initial conditions, then we reach the scaling regime of the attractor solution. However, if some of the fusion rates are reduced far enough to 'trap' a proportion of the DN in a small-k region for long enough, then we will not reach the scaling regime. Counter-intuitively, this can result in building up  Figure 4. Illustration of transition between scaling regime and bottlenecked regime. Left: Mass distribution at w = 25, starting from initial conditions of single nucleons, y k (0) = δ k1 , for the original kernel with K 1,1 = 4 (purple), and modified kernels with K 1,1 = 10 −4 (red), and K 1,1 = 10 −5 (yellow). We see that, within this range, making 1 + 1 fusions slower results in building up larger DN. Right: Number distribution at w = 25, again starting from initial conditions of single nucleons, for K 1,1 = (4, 1, 10 −1 , 10 −2 , . . . , 10 −6 ). The dotted line aligned with the K 1,1 = 10 −6 (green) curve is ∝ k −2/3 . We transition from converging very quickly to the scaling solution, to ending up with a power-law distribution cutting off at larger k. Between K 1,1 = 10 −5 and 10 −6 , the maximum k reached hardly increases, with the overall number density in the power-law tail just going down-we have reached the freeze-out limit for building up large DN.
larger DN than would otherwise have been the case. As roughly described by eq.(2.3), this occurs because small + large fusions are less velocity-suppressed in kinetic equilibrium than large + large fusions, so, if there is a bath of small DN present throughout the aggregation process, build-up interactions may freeze out at higher A. Figure 4 shows a simple example of moving between the scaling and bottlenecked regimes, which may be helpful to keep in mind through the following.
Since small + large fusions keep the number density of large DN the same, they act to increase the rate of large + large fusions, which goes as A 2/3 A −1/2 = A 1/6 . If Γ/H becomes of order 1 or higher, then large-large fusions will start operating and bring it down to ∼ 1, establishing a scaling distribution for the larger DN. Since the rate of fusions for a single small DN is larger than that for large DN, this also means that all of the small DN will be used up, clearing the bottleneck and placing us in the scaling regime.
The more interesting case is when large-large interactions are always frozen out. We can then model the aggregation process as a combination of the slow creation of large post-bottleneck 'seed' DN, and the fast accretion of small pre-bottleneck DN onto these. Looking first at the accretion process, from our previous assumptions the growth rate for a large DN scales geometrically, ∝ R 2 ∼ A 2/3 . So, dk/dw ∼ k 2/3 y s , where y s parameterises the concentration (and size) of small DN, and thus k ∼ dw y s 3 . If the bath of small DN is populated throughout most of the w time, then k max ∼ w 3 max , realising the freeze-out bound of eq.(2.3).
As well as the maximum size, we are also interested in the number distribution over k. If a seed DN is produced at time w inj , then its eventual size will (for roughly constant y s ) be k ∼ (w max − w inj ) 3 . More precisely, if we change time variable to z, where dz/dw = y s , we have k ∼ (z max −z inj ) 3 . This means that the relationship between injection time z inj and final size k is given by − dz inj dk ∼ k −2/3 . Next, define the injection rate f i ≡ dy inj dz , where y inj parameterises the concentration of injected seed DN. Since accretion only changes the knumber, and not the density, of an injected seed population, we obtain dy k dk ∼ k −2/3 f i (z inj ). In particular, if for small z the injection rate is roughly constant (this is reasonable, since only a small fraction of the pre-bottleneck bath has been used up, and the temperature, scale factor etc. change by only a small amount), then for large k we should have a power law number distribution dy k dk ∼ k −2/3 . As we increase the injection rate f i , we move through the three different regimes described at the end of Section 2.1: • If the injection rate is small enough, then most of the small bath is not used, and most of the mass is stuck behind the bottleneck, with a small proportion in a dy k /dk ∼ k −2/3 f i (z inj ) tail, which extends to k max ∼ w 3 max . As an example, if we assume that f i decreases with time, such that the nucleon number integral kmax dk k k −2/3 f i is dominated by the upper decade, ∼ k 4/3 max f i (0), then the upper limit of this regime is given by • For larger injection rates, we use up all of the small bath before we reach w max , so we build up to correspondingly smaller k max . Under the same assumptions on • For sufficiently large injection rates, large-large fusions eventually become important, and we enter the scaling regime. We expect this to happen roughly when k max for the addition process is comparable to the k values of the scaling peak. Since this has k ∼ w 6/5 , under the above assumptions on f i we expect the cross-over to be at In the first regime, the proportion of DM in the pre-bottleneck bath stays roughly constant, so dz dw = y s const. From the previous section, ∆w ∼ ∆(a −(1+γ) ), so the bulk of the process again takes of order a Hubble time. For the intermediate regime, we use up most of the small bath within ∆w (3k max ) 1/3 < w max , a small fraction of a Hubble time. For the rest of w time, either large-large fusions are frozen out, in which case we only make small modifications to the number distribution, or we enter the scaling regime. Figure 4 shows the numerical solution of a particularly simple bottlenecked exampleusing the geometrical kernel of eq.(3.5), but reducing K 1,1 . In this case, the injection rate f i ∝ y 2 1 , so for most of the addition process it is constant, giving a power-law number distribution dy k /dk ∼ k −2/3 for large k. The right panel of Figure 4 illustrates how we move through the three regimes identified above as we decrease K 1,1 . We start out converging quickly to the scaling distribution, but for K 1,1 10 −4 we never reach this regime, ending up with a power law out to larger k. For K 1,1 10 −5 , most of the k = 1 population never makes it past the bottleneck, and our power-law tail goes to smaller concentrations rather than higher k values. The right panel of Figure 5 shows explicitly that, for the bottlenecked solutions, we are close to being in the addition-dominated regime, i.e. the dominant process is 1+k → (k+1). Comparing the left panel to Figure 4 illustrates that, as we increase w max , the difference between the k max attainable in the scaling and bottlenecked regimes increases.

Aspects of dark sector phenomenology
So far, we have investigated the BBDN process, and the number distribution of DN that may result from it. The DM self-interactions, and the extended hidden sector needed to realise such models, mean that these theories generically have the possibility of interesting hidden sector phenomenology, including non-standard indirect detection signals, modifications of halo properties, and early-universe signatures of extra species. Also, though not required in these models, it is possible that the DN have sufficiently strong interactions with SM states to give signals in future direct detection experiments. If that is the case, then such signals may differ from those of usual WIMP DM, and have different relationships to collider bounds [53]. Additionally, such interactions may lead to the capture of DN in astrophysical objects, and then self-interactions among captured DN could become important.
In standard ADM models, annihilations with the relic symmetric population can lead to astrophysical energy injections, most detectably in the early universe, putting a limit on how large this population can be. In nuclear DM models there is additionally the energy injected by inelastic processes, particularly the release of binding energy from fusions.
For ADM of mass 100 GeV, annihilation of the symmetric component generally needs to go into, or via, lighter hidden sector particles rather than directly to the SM, to avoid constraints from direct detection and collider experiments [48]. We will see below that in the case of DN, the constraints from direct detection on the SM couplings of individual DN constituents are often even tighter, confirming the need to annihilate into lighter hidden sector particles. In addition, if we want to build up DN via fusions from single nucleons, we need some lighter state for at least the first few fusions to de-excite into, as the limits on SM interactions are generally strong enough that we cannot de-excite fast enough via those alone. We also need some mediator particle to transmit the binding force. If any of these additional states are sufficiently light, strongly-interacting, or long-lived, then there may be astrophysical constraints on their properties.
The amplitudes for sufficiently weak interactions between DN and other states will, for interactions involving individual dark nucleons, factorise into a form factor and a singlenucleon amplitude. For momentum exchanges which are small compared to the inverse radius of the DN, the form factor will be ∝ A, i.e. the interaction will be coherence enhanced, and will fall off in some nuclear-structure-dependent way for larger momenta [53] analogously to the scattering of SM nuclei with WIMP DM.
Assuming that DN scatterings with SM nuclei in DM direct experiments are in the small momentum transfer regime, the coherence factor has the effect of enhancing, for a given SM-dark nucleon interaction strength, the constraints from these experiments. Taking the DN to be of a single size A, the number density is ∝ 1/A, but the scattering rate is coherence-enhanced by A 2 , so the overall event rate scales as A. For the distributions discussed in the previous sections, both the mass distribution ky k and the scattering-rate distribution k 2 y k are peaked over a small logarithmic range in k, so the ∝ A enhancement is a reasonable approximation. The same is true, for simple injection profiles, in the bottlenecked case. Since the rate for processes with only SM particles in the initial state is generally set by the interaction strength with individual dark nucleons, collider constraints, and those from the annihilation of the symmetric DM component in the early universe, are relatively less stringent compared to direct detection.
In addition to these effects on overall rates, for intermediate momentum transfers the effect of the momentum-dependent DN form factor may be important. The possibility of 'dark form factors' influencing the shape of direct detection recoil spectrum has been investigated by a number of authors, often in the context of suppressing low-momentum scatterings [60,61], but also addressing soft form factors [62,63]. Also, for sufficiently large energy transfers, inelastic processes involving excited DN states or fissions may be important. In a forthcoming paper [53] we investigate the possible direct detection phenomenology of DN in quantitative detail, as well as exploring the possible consequences of its capture in astrophysical objects.

Post-nucleosynthesis energetics in the dark sector
Since we have been considering DM with strong self-interactions, an obvious question is whether these self-interactions have late-time consequences. For elastically-scattering DM with velocity-independent self-scattering cross section, the observational limit on the selfscattering cross section is σ XX /m X 1 barn/ GeV (see e.g. [64]). For a population of in the notation of eq.(2.2). From the results of Section 2.1, if we scale M 1 → λM 1 , ρ b → λ 4 ρ b etc., then in the scaling regime, A → λ −12/5 A. Thus, σ AA m A → λ −11/5 σ AA m A , so we would hit the observational limit at ∼ 10 MeV scales for the constituents, assuming that we built up to the largest allowed sizes.
Realistically, the situation can be more complicated than this. Firstly, there would a spectrum of DNN values. However, for the scaling distribution, and simple forms of the bottlenecked distributions, the area distribution k 2/3 n k and the mass distribution kn k both receive most of their contribution from a single decade or so of nucleon numbers, so the above calculations will hold approximately.
More interestingly, the fact that DN collisions may not be elastic could have phenomenological consequences. Since, neglecting collisions, DN of different sizes will have the same velocity distribution, with e.g. v ∼ 10 −3 in the Milky Way, large DN will have kinetic energies much larger than the single-nucleon binding energy scale. Thus, it is possible that DN-DN collisions will result in inelastic collisions, including nuclear fragmentations, transitions to excited states, and fusions.
Some types of collisions may be 'endothermic', reducing the average kinetic energy (KE) per nucleon of the final state DN. Examples include scattering to excited states, in which some of the initial KE is lost into de-excitation products (fusions with sufficiently small binding energy differences also lose KE overall). If these de-excitation products interact with the SM, there may be indirect detection signals from these processes. Independently, the inelastic form of the scattering will tend to lead to the contraction of DM mass distributions, as compared to elastic scatterings. It may be the case that, in regions of sufficiently high DM density, there is the possibility of run-away contraction (in particular, if the DM is self-gravitating, then removing KE results in the distribution contracting and also heating up, due to the negative heat capacity). More generally, the effects and constraints will be different from the usual simulations of elastic DM self-interactions.
On the other hand, the presumed binding-energy-based stability of large DN means that there may be exothermic collisions, in which the average KE of the DN is increased (c.f. fusions in stars, in the SM). If the velocity kick given to the DN is large compared to the velocity dispersion in the halo, exothermic collisions will have the effect of clearing out the central high-density region, until the number density is low enough that collisions occur once per particle per Hubble time. Since smaller halos tend to have smaller velocity dispersions, small velocity kicks may only modify structure on small scales (where there are problems with Λ-CDM predictions for structure).
These general effects of inelastic DM self-interactions are not specific to nuclear DM, and a number of models have been investigated in the literature. In particular, models in which the DM transforms under some approximate symmetry can easily realise small mass splitting, and the effects of illustrative cases corresponding to a Yukawa self-interaction have been considered in [65,66].
Turning to the energy released in de-excitations; considering energy injections once the aggregation process has frozen out, the injection rate is just set by the local DM velocities and number density (which, in regions where only a small fraction of the DN undergo self-interactions, will be the standard collisionless DM values). Generally, if the rate is velocity-suppressed (as for large DN-DN collisions), and freeze-out is before BBN, then present-day cosmic ray (CR) constraints dominate those from earlier times, if the injection is to sufficiently energetic SM particles [67,68]. Velocity-independent processes with rate ∝ n 2 have BBN, CMB and CR constraints within a few orders of magnitude of each other, depending on the form of SM injections [67][68][69].
In regions where only a small fraction of DN undergo fusions, the proportion of the DM mass density represented by the binding energy released in late-time fusions is where we take t gal ∼ 10 Gyr. Since the detectability of CR signals from annihilating DM at ∼ pb cross-sections depends on the SM injection channel etc. [70], the same applies for DN collisions within the A ranges of interest. Given the range of initial DN masses, and the possibility of excited states of DN, such signals may have a richer structure than the indirect detection signals usually considered. Additionally, the geometric cross sections we have assumed give different velocity dependence to the partial wave processes usually considered. Also, as considered above, for endothermic collisions in sufficiently dense regions, there may be the possibility of (run-away) contraction of the DN distribution, which could significantly affect signals from e.g. the Galactic centre.

Light dark sector states
As mentioned above, DN models generally require additional lighter hidden sector states. Since the assumption is that these do not make up the bulk of the DM, they either need to never have a large abundance (generally difficult to realise in thermal ADM histories), be sufficiently light and weakly interacting with the SM to persist as dark radiation, or their abundance needs to be reduced.
Reducing the yield of a species can occur by transferring its energy density to other hidden sector species, or to the SM. In the latter case, if this injection occurs during or after BBN, there are constraints on its form and magnitude. For injection to hadronic channels, the total energy density injected after T γ ∼ 1 MeV must be significantly below that of the DM. For non-hadronic injection, the dominant effect at times before thermalisation becomes inefficient, 10 4 sec, is alteration of the photon-baryon ratio 11 , which constrains the amount of energy injected to 0.1ρ γ [71]. At later times, any energy injections, apart from those to neutrinos, must again be significantly sub-DM [72][73][74].
For hidden sector particles of mass 100 MeV, their chemical equilibrium abundance (at the SM temperature) by BBN times is significantly sub-DM. The dominant constraints on their couplings to the SM generally come from collider experiments, and permit decay times of 1 sec (e.g. [75]). So, in these cases, if such decays are possible, then there are generally consistent scenarios in which any initial abundance decays to the SM before BBN. 12 On the other hand, for hidden sector particles with m 100 MeV, astrophysical observations place strong constraints on their SM couplings. These restrict their SM decay times to 1 sec 100 MeV m 1+2k , where k ≥ 0 is set by the mass dimension of the decay operator (e.g. [76,77]), and mean that direct interactions with the SM are frozen out below very high temperatures (T γ m). If m 10 keV, purely SM decays occur at T γ < m, so the decay of a thermal abundance would transfer ρ γ energy density to the SM. For smaller m, the decay time is 10 4 sec, so the constraints on energy injection to the SM are much more severe. In either case, limits on SM energy injection generally imply that the majority of the initial thermal energy density needs to be transferred to other, lighter hidden sector states. 13 Considering these states in turn, we eventually require that there be long-lived hidden sector states, which will act as dark radiation, at least during the early universe. Extra relativistic species are compatible with current observations [78], and if the hidden sector is at a lower temperature than the SM, its contribution to the effective number of such species N eff may be small. Plausible candidates for these species within hidden sector models include very light pseudo-Nambu-Goldstone bosons or Z states. 11 and also increasing the photon temperature relative to that of neutrinos, so decreasing N eff at later times compared to BBN. 12 As noted in previous sections, obtaining sufficiently fast annihilations (directly to the SM) to reduce the energy density to sub-DM levels is generally difficult for m 100 GeV. 13 It may be possible to realise alternative scenarios in which energy is transferred to the SM indirectly, via other hidden sector states. If the hidden sector is at a slightly lower temperature than the SM, and the transfer is before ∼ 10 4 sec, this may be safe (sitting at the upper end of the mass range, and transferring mostly before BBN, may also work). This would require that number-changing interactions with other hidden sector states are sufficiently fast to keep the light state in thermal equilibrium during the process. Three-body interactions with heavier hidden sector states are suppressed by two factors of the heavy state number density (which is small for DM of SM mass or higher), and inelastic two-body collisions are only relevant if such excitations are accessible at low temperatures, imposing model-building constraints.

Conclusions
In this paper we have studied the "Big Bang dark nucleosynthesis" process by which 'nuclear' bound states of DM may be built up in the early universe. Specifically we focussed on the case of asymmetric DM models where the nuclear binding energy per dark nucleon saturates in the large nucleon number limit. We find that, if fusions between small dark nuclei (DN) happen sufficiently fast, and fusion cross-sections between large DN scale on average geometrically, the resulting number distribution generically takes on a universal form, illustrated in Figures 1 and 2. This result is broadly independent of small changes to the initial conditions or fusion rates, assuming that the DN are in kinetic equilibrium throughout the process, as discussed in Section 3.2. The average mass of the DN built up during this process can be as large as ∼ 10 8 GeV.
If fusion reactions between small DN are not large enough to reach this regime, but large-small fusions still occur at an appreciable rate, then there is the counter-intuitive possibility of building up even larger DN due to the higher velocities of smaller particles, leading to larger fusion rates. The resulting number distribution generally takes the form of a power law multiplied by an 'injection profile' parameterising the change of small-small fusion rates with time. Figures 4 and 5 illustrate this behaviour.
Given this possibility of building up large dark nuclear bound states, direct detection signals may be modified in a number of ways-by the coherent enhancement of SM-DM interactions, by dark form factors if the DN radius is significantly larger than SM nuclear radii, and by inelastic interactions if there are sufficiently low-lying excitations. We reserve detailed discussion of such direct detection phenomenology to a forthcoming paper [53].
The possibility of inelastic collisions can also lead to interesting astrophysical dark sector interactions as mentioned in Section 4.1. These include annihilation-type indirect detection signals, which are usually absent from asymmetric DM models (and here may have a richer structure, corresponding to the number distribution of DN). Inelastic collisions may also modify the effects of DM self-interactions on halo structure, especially at short distance scales. In addition, as discussed in [53], there are a range of possible consequences for the capture of DN by astrophysical objects as well, from ejection by de-excitations, through to fusions leading to a very dense DN core.
Finally, the discussion in this paper has intentionally been as model-independent as is practical, investigating idealised versions of the behaviour that classes of models can display. The physics is inspired by that observed in the Standard Model (SM) sector, with the simple, but important, change that the Coulomb barrier is removed. Although the example of the SM reassures us that a sufficiently complicated model can realise the physics we discuss, it would be interesting to investigate, along the lines of [9,14,79], specific simple toy models which realise all or some of the features discussed in this paper. As with the SM, we would expect there to be additional features in the dark nuclear spectroscopy and interaction cross sections, e.g. due to shell structure, on top of the general scaling behaviour of nuclear properties that we have utilised. In addition the extra hidden sector states required (to mediate binding forces, carry away de-excitation energy, etc) beyond the dark nucleon itself can be relevant for astrophysical phenomenology and constraints.
In particular, if some of these states have masses 100 MeV, there is likely to be a need for extra hidden sector states which act as dark radiation. Another avenue for further work is the possibility of symmetric DM models.

A Transition from equilibrium to aggregation
In chemical equilibrium, detailed balance for the k + (A − k) ↔ A process implies that Γ A→(k,A−k)ñA = σv (k,A−k)→AñkñA−k , (A.1) where theñ k are the equilibrium number densities, Γ A→(k,A−k) is the rate of A → (k, A−k) dissociations for a single A-DN (in general, this depends on the abundances and velocities of other species), and σv (k,A−k)→A is the thermally average fusion cross section (we assume, as is almost always the case, that mean free paths are long enough so that velocity rather than diffusivity is important). If we assume that the fission rate is purely determined by the temperature, then when we are in kinetic equilibrium at (non-equilibrium) number densities n k etc., This assumption about the fission rate will be violated if DN-DN collisions often lead to 'prompt' fragmentations, rather than either elastic scatterings, or inelastic collisions that result in a thermally-equilibrated excited state (c.f. the 'compound nucleus' model for SM nuclear collisions [54]). The same is true for DN collisions with other baths that do not have their chemical equilibrium abundance. Assuming ideal gas behaviour, the equilibrium number densities are related bỹ where g k is the number of effective degrees of freedom of k-DN. Writing m k = km 1 − B k , and assuming that binding energies are a small fraction of the mass, this is where Λ = 2π/(m 1 T ).
The overall forward rate for the k + (A − k) ↔ A process is σv n k n A−k − Γn A = σv n k n A−k 1 −ñ kñA−k n A n A n k n A−k . (A.5) Thus, as long as dissociations will give only a small correction to the fusion rate. For the 1 + 1 ↔ 2 process, we start building up when n 1 Λ 3 e B 2 /T > 1. More generally, if we assume that typical binding energy differences rise fast enough with size difference (for example, in a liquid-drop type model B k = αk − βk 2/3 , the surface tension term β gives B A − B A−k − B k ∼ βk 2/3 for A > k), then if n 0 Λ 3 e α/T 1 (where n 0 is the number density of nucleons), the inequality should hold for all of our fusion processes between k, A − k states with reasonable density, with the possible exception of anomalously low (or wrong-sign) binding energy differences for small-number states.
Since n 0 Λ 3 1 (the DM is a dilute gas), we have equality n 0 Λ 3 e α/T = 1 at T α. After that, assuming that −d log T d log a = O(1), we have d log LHS d log a ∼ α T 1. So, with the possible exception of small small-k B k values, we go from being near-equilibrium, to dissociations being very suppressed, in a time that is parametrically H −1 T α . Furthermore, during this time the build-up of DN proceeds at most as quickly as a fusions-only aggregation process, like that considered in Section 3. As explained there, in many regimes this progresses fairly uniformly over the first Hubble time, so a modification to the fusion rates during a small initial fraction of this time will have a small overall effect on the result.