Early Universe synthesis of asymmetric dark matter nuggets

We compute the mass function of bound states of asymmetric dark matter — nuggets — synthesized in the early Universe. We apply our results for the nugget density and binding energy computed from a nuclear model to obtain analytic estimates of the typical nugget size exiting synthesis. We numerically solve the Boltzmann equation for synthesis including two-to-two fusion reactions, estimating the impact of bottlenecks on the mass function exiting synthesis. These results provide the basis for studying the late Universe cosmology of nuggets in a future companion paper.


I. INTRODUCTION
Searches for dark matter (DM)-whether through DM annihilations to Standard Model (SM) particles in our Galaxy, through observations of structure, through collider searches, or through DM interactions with SM nuclei in underground detectors-have focused on DM that behaves as a single massive particle or as a coherent field configuration. The weakly interacting massive particle (WIMP) and the axion have served as primary motivators of these search techniques; in addition to being well-motivated candidates, they provide sharp predictions for experiments probing the nature of dark matter. It is, however, important to explore well-motivated new ideas for dark matter candidates that lead to qualitatively different experimental signatures of dark matter, and inform new search strategies.
Even modest changes in the nature of the dark sector can have radical implications for cosmology and its associated observational constraints. In this paper, we consider a DM particle with a particle-antiparticle asymmetry [as in asymmetric dark matter (ADM); see e.g. [1,2] for reviews] that self-interacts through an attractive force. Similar to SM nuclei, given a sufficiently strong attractive force, DM particles bind together to form composite states. The size of these composite states depends both on the strength and range of the force, and on the presence or absence of bottlenecks. For example, in the SM, only small nuclei are synthesized in the early Universe because the A ¼ 5, 8 nuclei are unstable; three helium nuclei must fuse to carbon in order to surpass this bottleneck. In a dark sector, however, such bottlenecks may not be present. This is especially true if a repulsive long-range force like electromagnetism is absent, and if the force mediator binding the dark nuclei is longer range than the nuclear force mediated by QCD mesons.
We study the synthesis of the dark nuclei, which, as in Ref. [3], we refer to as "dark nuggets (DN)." Our goal is to obtain results that connect the size of the synthesized nuggets to a simple UV complete model; we employ a model with fermionic ADM and a scalar mediator with arbitrary mass that mediates an attractive dark force. We draw on the results of our companion paper [4], which computes the relevant nugget properties-notably the saturation density and binding energy-using the σ-ω model from nuclear physics. The σ-ω model, named for the lightest scalar and vector QCD mesons mediating the most important attractive and repulsive interactions that bind large nuclei, employs relativistic mean field theory (RMFT) in order to solve the equations of motion that determine bulk properties of ground states of nuclei. We use results from [4] along with the compound nucleus model to obtain approximations for the fusion and splitting cross sections.
Synthesis becomes efficient once the temperature of the Universe drops below the binding energy of the nugget. Utilizing simple analytic arguments, along with our results in Ref. [4], we are able to obtain an estimate for the typical nugget size just before the era of structure formation. We also obtain a distribution of nugget sizes by solving the Boltzmann equations for synthesis numerically. Using the results of this simulation and argument by analogy to SM nucleosynthesis, in our Coulomb repulsion-free scenario, we argue that a substantial bottleneck is likely to occur only if multiple states adjacent in size (e.g. N ¼ 3, 4) are unstable. Superficial bottlenecks, corresponding to isolated states being hard to form (either due to a very low formation cross section or absence of a stable state), do not significantly affect the distribution of nugget size (or mass function) at the end of synthesis.
Our main results are shown in Figs. 2 and 5, for the cases of a presence or absence of a bottleneck to synthesis at low N. In parameter ranges where two-particle bound states are easily formed in the early Universe (see [3]), assuming that other substantial bottlenecks at low N do not occur, large nuggets can be synthesized in the early Universe in the presence of sufficiently attractive and long-range forces. For example, as shown in Fig. 2, assuming DM mass m X ∼ GeV, DM-mediator fine structure coupling constant α ϕ ∼ 0.3, and mediator mass m ϕ ∼ MeV, the typical nugget size exiting early Universe synthesis is of order N ∼ 10 12 or M N ∼ 10 11 GeV, with the sizes scaling as N ∼ ðBE 3=2 2 =m 7=2 X Þ 6=5 and M N ∼m X N, where the twoparticle bound state binding energy that sets the synthesis temperature is BE 2 ∼ α 2 ϕ m X =4 and the nugget saturation mass scale ism X ∼ α −1=4 ϕ ffiffiffiffiffiffiffiffiffiffiffiffi ffi m X m ϕ p . Since large nugget binding energies are typically quite large when the DM/ mediator mass hierarchy is large, the nugget masses can be significantly smaller than m X N.
If synthesis proceeds at all, we find that nuggets are generically synthesized past their saturation size, where the force range, m −1 ϕ , sets the nugget size and binding energy along with the coupling and DM mass. Nuggets bound by a scalar mediator much lighter than the dark matter constituent saturate at a larger N, but by virtue of the effectively stronger coupling in this limit, synthesis to large nuggets occurs more readily. Simulating synthesis with fusion through both mediator and small nugget emission, we show that the spectrum of nugget sizes exiting synthesis can be fairly broad. The mass function exiting synthesis is shown in Fig. 4. This may have implications for the late Universe cosmology and detection prospects, which we study in a future companion paper [5]. Substantial bottlenecks at low N, on the other hand, can lead to a bimodal distribution of dark matter exiting early Universe synthesis, with the majority of DM residing in the form of small-N nuggets and a subdominant population of even larger nuggets-for example, as shown in Fig. 5, with typical size N ∼ 10 20 or M N ∼ 10 19 GeV for the benchmark quoted above.
Dark matter nugget properties and synthesis ("darkleosynthesis") have been explored elsewhere, though with more limited assumptions for obtaining quantitative results. Early Universe synthesis of two-particle bound states was addressed in [3]; we will utilize these results in the initial step of our analysis. Reference [6] considered larger N bound states, but did not address the saturation limit. We will find, however, that nugget synthesis typically terminates with saturated nuggets, even when the attractive force is very light relative to the DM. Reference [7] considered the synthesis of nucleuslike nuggets, assuming dark nuclei properties (notably the saturation density and binding energies) that directly mirror SM nuclei. We utilize similar analysis techniques and confirm several major results of that work-that the typical size of synthesized nuggets scales as a dimensionless interaction time to the 6=5th or 3rd power in a no-bottleneck or bottlenecked scenario, respectively, and that, in the no-bottleneck scenario, nugget sizes are peaked within a couple of generations in size about the typical size. Our work differs from that of [7] in two primary ways. First and most importantly, we provide a detailed accounting of underlying assumptions entering the setup of Boltzmann's equations within a concrete model through employing the results of [4] for nugget properties, which also allows us to link the dimensionful scales entering the interaction time to Lagrangian parameters; this reveals that larger synthesized nugget sizes (by several orders of magnitude) than advertised in [7] are possible. At the same time, it is notable that our analysis justifies many of the arguments made in [7]. Second, instead of employing a simple fusion model with geometric cross sections involving a single mediator emission, we consider the compound nucleus model and show that multiple mediator emissions are dominant, but also consider the effects of possible light dark nugget emission, analogous to nuclear fusion through nucleon or α particle emission. Reference [8] studied the synthesis of dark nuclei modeled on SM nuclei, with a dark confining force binding the composite dark baryons into nuclei and an additional weak dark electromagnetism enabling "di-darkleon" fusion. Lastly, Ref. [9] studied the synthesis of dark spin-0 deuterium forming in a two-flavor, two-color, dark QCD.
The outline of this paper is as follows. In the next section we review the features relevant for synthesis of our nuclear model for ADM nuggets. Then in Sec. III we outline the conditions for synthesis to begin in the context of our simple UV complete model. In Sec. IV, we set up the Boltzmann equations with the appropriate fusion and dissociation rates before solving them. We utilize these results to obtain analytic estimates for the typical nugget size. We conclude with an eye toward future work exploring the impact of dark nuggets on stellar and structure formation.

II. A NUCLEAR MODEL FOR ASYMMETRIC DARK NUGGETS
We will consider a model with a single Dirac fermion with attractive self-interactions mediated by a lighter, real scalar, governed by As discussed in detail in [4], large collections of dark matter can form stable bound states when α ϕ . Two-body bound states form when α ϕ 2 m X m ϕ ≳ 0.84 (assuming perturbative coupling; see [3] and references therein). Here we are interested in the synthesis of such bound states in the early Universe. As we will see, for synthesis to begin with two-particle bound states, it is important that the force mediator be light enough to be produced on shell in the fusion process, such that BE 2 ≳ m ϕ . Since, for perturbative coupling, BE 2 ∼ α 2 ϕ m X , one generally needs m X > m ϕ .
It is also natural to consider a strongly coupled dual of this model, where X is the analog of a nucleon and ϕ is the analog of the lightest scalar meson, σ (or f 0 ). In a composite model, one expects other meson degrees of freedom such as vectors (like the ω) and pseudoscalars (like pions) to mediate repulsive and/or spin-dependent interactions of comparable importance. If there is an additional approximate symmetry in the dark sector, the analog of isospin-dependent interactions may also be important. One also naively expects a mass hierarchy between the dark matter constituents (baryons) and force mediators (mesons) to be moderate, and for the masses of the mesons to be of very similar size. As discussed in [4], if this is the case, it is unlikely that there is a viable region of parameter space in which any of the mesons is lighter than the two-body binding energy, which will stifle early Universe synthesis.
In the SM, deuterium forms through p þ n → D þ γ. The deuterium binding energy ∼2 MeV is significantly smaller than m π 0 ∼ 135 MeV. Without electromagnetism, fusion into deuterium would require π 0 emission. Such a process would then only be efficient at temperatures near m π , where dissociation would dominate. So, taking a cue from the SM, why not just add a dark photon in order to enable the first step of synthesis? (This is precisely the scenario considered in [8].) There are costs. Including a dark photon will destabilize dark nuggets of sufficiently large size-just as electromagnetism helps to destabilize large nuclei. The larger the coupling, the smaller the size at which nuggets will destabilize. So allowing for very large nuggets requires very small coupling. But the smaller the coupling, the smaller the two-body fusion cross section.
Due to these complications, we will focus on the case of an attractive force mediator only, where there is a large parameter space for efficient nugget synthesis. We will restrict our attention to perturbative couplings, where we have a good handle on the two-body bound state properties that govern the initiation of early Universe synthesis, even though the RMFT used to deduce properties of large nuggets is valid also for nonperturbative coupling. In our synthesis estimates that follow, the behavior of binding energy and fusion cross sections as a function of nugget size will be important. Thus here we first summarize these features of nuggets, justified and presented in more detail in [4].

A. Binding energy and the liquid drop model
In [4], relativistic mean field theory was used to derive the behavior of nugget structure (density, size) and binding energy as a function of the nugget number, N. RMFT applied to nucleons has been used to accurately model bulk properties of large nuclei such as binding energy, density, and the saturation property of nuclei-that the density is relatively constant as a function of mass number [10].
As one expects since the only force in play is attractive, the binding energy per dark number increases as a function of N. At some N determined roughly by N sat ∼ ðα ϕ , the binding energy per dark number levels off, asymptoting towards a constant determined as a fraction of m X by the combination of parameters, C 2 ≡ α ϕ m 2 X m 2 ϕ (a larger fraction for larger C 2 ). At N > N sat , the density as a function of N also approaches a constant, exhibiting saturation behavior: adding further constituents does not change the nugget density but simply increases the nugget size as RðNÞ ∝ N 1=3 . For N > N sat , just as for large nuclei which exhibit saturation, a liquid drop model gives a good description of the nugget binding energy: where 0 < μ 0 < m X is a bulk energy constant equivalent to the chemical potential of infinitely large nuggets and is a correction term proportional to the surface area of the nugget that takes into account the dearth of close-range attractive interactions of constituents at the surface which would have decreased the energy of the configuration. The infinite matter chemical potential μ 0 approaches m X as C 2 → ∞ and the surface energy constant ϵ surf grows with decreasing m ϕ =m X and decreasing α ϕ . (See Table I of [4].)

B. Saturation properties of nuggets
In the following sections, we will show that in the range of parameter space in which two-body synthesis proceeds easily, nuggets are generically synthesized up to sizes beyond which they exhibit saturation behavior. We are able to self-consistently estimate the results of synthesis in terms of the saturation properties of nuggets: chiefly, in terms of average nugget mass per dark number,m X ¼ m X − BE N =N ∼ μ 0 and nugget length scale r 0 ≈ ð 4 3 πn sat Þ −1=3 , with n sat the saturation number density.
In the limit where large nuggets are strongly bound (which is generally true when α ϕ < 1 and BE 2 > m ϕ ), the saturated nugget parameters r 0 andm X0 are given approximately bȳ Them X0 → m X limit represents the limit of no binding, where the above equations are invalid. Inclusion of a quartic potential of the form VðϕÞ ¼ 4 3 λ 4 α 2 ϕ ϕ 4 leads to an effective scalar mass, where m Ã is the effective DM mass, related to the scalar field VEV through m Ã ¼ m X − g ϕ hϕi. The equations of motion guarantee that m Ã ≥ 0. In the limit where α ϕ m 2 X ≫ m 2 ϕ and λ 4 ≪ 1, m Ã ≪ m X and we find Eq. (3) holds if we replace m ϕ with m ϕeff . See [4] for details.
The saturation size, N sat ∼ ðr 0 m ϕeff Þ −3 , corresponding to the point where nugget radius exceeds the force range, R ≳ m −1 ϕeff , is also important in checking the consistency of our estimates. In the synthesis estimates that follow, we will assume λ 4 ≪ 1 so that the effects of the quartic potential are essentially encapsulated in m ϕeff . We will frame our analysis as if the potential is completely absent, but we note that our results hold equally well when a quadratic potential with λ 4 ≲ 0.01 is included through taking m ϕ → m ϕeff . It should also be noted that, even if a large hierarchy between m ϕ and m X is achieved, the quartic coupling will limit the size of the m X =m ϕeff hierarchy relevant to nugget properties.

III. CONDITIONS FOR SYNTHESIS
Here we discuss conditions for initiating synthesis with the formation of two-dark-nugget bound states, 2 X, the temperature at which this initiation occurs, and possible bottlenecks at low dark nugget number. We will follow nuclear physics convention and denote each dark nugget species as N X, where N is the dark nugget number.

A. Conditions for beginning synthesis
Nugget synthesis begins in the early Universe when 2 X bound states form. This very first stage of synthesis corresponds to passing the two-DN bottleneck, analogous to the deuterium bottleneck in big bang nucleosynthesis (BBN). The bound state 2 X can begin to accumulate if the 2 X formation rate exceeds the Hubble expansion rate when the dissociation rate drops below the formation rate. This occurs when the number density of mediators energetic enough to dissociate the 2 X drops below the 2 X equilibrium number density.
Because formation through the process X þ X → 2 X þ ϕ is typically inefficient when m ϕ ≳ BE 2 , as argued in detail in Ref. [3], efficient two-DN synthesis generally requires where we have used the expression for the two-body Yukawa bound state binding energy in the Coulomb (hydrogenlike) limit, which is self-consistent as long as the coupling is perturbative. 1 In addition, for the formation rate to exceed the Hubble rate at any temperature below the X freeze-out temperature, we also require (see [3]) ðsecond 2 X synthesis conditionÞ: ð6Þ Figure 1 shows these 2 X synthesis conditions, along with the region of parameter space for which the 2 X synthesis temperature-as we discuss in the next section-is an order of magnitude larger than the temperature at matter-radiation equality, T syn > 10 T eq . Synthesis near matter-radiation equality would require an analysis taking into account nonlinear structure formation and is thus neglected. It is interesting to note that synthesis can begin significantly before the end of radiation domination only if α ϕ ≳ 0.01.

B. Synthesis temperature
Since dissociation of 2 X must become inefficient relative to formation in order for synthesis to begin, and because dissociation becomes inefficient only once the temperature drops below the binding energy, the big bang darkleosynthesis temperature, T syn , is typically set by BE 2 . This is analogous to BBN, where T BBN is set by the deuteron binding energy, BE D ∼ 2 MeV.
FIG. 1. Approximate conditions for synthesis to begin with efficient formation of 2 X. The top line corresponds to the condition for the formation rate to ever be larger than the Hubble rate in the BE 2 ≫ m ϕ limit [Eq. (6)]. The bottom line shows the condition for synthesis to begin when the temperature is an order of magnitude larger than T eq , ensuring that our analysis assuming synthesis before the onset of nonlinear structure formation remains valid. The dashed lines correspond to BE 2 ¼ m ϕ ; at least roughly speaking, BE 2 > m ϕ [Eq. (5)] in order for the dissociation rate to ever dip below the formation rate. 1 Thermal corrections to m ϕ may be relevant, but since they depend on interactions that do not directly contribute to fusion, we assume that they are small and can be neglected.
More specifically, one can estimate T BBN by setting the number density of photons with energy greater than BE D equal to the deuteron equilibrium number density, and solving for the temperature. Because the baryon-to-photon ratio is so small, the temperature must fall to order T BBN ∼ BE D =37. Following the same logic, and ignoring all but the 2 X ground state, we see that chemical equilibrium implies a Saha relation, where we have taken the ϕ chemical potential to be zero, which is valid as long as number changing processes for ϕ are still efficient. At the very beginning of synthesis the total dark number density n X is dominated by unbound X's so that we may take n1 X ≈ n X ¼ Ω DM Ω B m p m X n B wherem X is the average mass per dark number at the end of synthesis.
The equilibrium number density of ϕ's with enough energy to dissociate 2 X is where again we have set μ ϕ ¼ 0 and T ϕ ¼ T X . Setting n ϕ ðϵ > BE 2 Þ ¼ n2 X leads to the condition for synthesis temperature, where we have assumed BE 2 ≪ m X . For MeV ≲ m X ≲ TeV, 0.01 ≲ α ϕ ≲ 0.3, and including other constraints for 2 X synthesis discussed above [Eqs. (5) and (6)], we find that BE 2 =T syn X ∼ 10 to 30. Here we briefly discuss the relationship between the DM and SM bath temperatures. Suppose the dark sector kinetically decouples from the SM at temperature T d , when the DM bath is still relativistic. Then, while the DM bath is still relativistic, the ratio of dark to SM temperatures changes only when some relativistic species falls out of chemical equilibrium and dumps its entropy in one bath or the other. Supposing only the SM bath is heated relative to the DM bath, for example, then counts the number of relativistic degrees of freedom in the SM when the two sectors decouple and g SM ðT γ Þ counts them at temperature T γ . When ϕ becomes nonrelativistic and decays, its entropy is dumped either into some hidden sector bath (which would heat that bath relative to the SM bath) or into the SM bath, heating it further. If m ϕ ≳ 5 MeV when it decays to the SM (whether to electrons or neutrinos), there is little effect on BBN or cosmic microwave background (CMB) measurements of the number of relativistic degrees of freedom, N eff [11]; if, on the other hand, m ϕ ≲ 5 MeV and the decay is to neutrinos or to additional radiation in the dark sector, there may be observable signals in next-generation CMB experiments, depending on the decoupling temperature of the dark sector from the SM. In this paper we take an agnostic view towards the exact dynamics that determines the ratio of temperatures, but note that in many plausible scenarios, T γ T X ∼ few, and we take T γ T X ≈ 1 in several illustrative plots that follow; this ratio has little quantitative or qualitative impact on our results.
Once T ≲ T syn X , larger N nuggets will efficiently form unless there is an additional bottleneck at low N analogous to the A ∼ 5 bottleneck for BBN. This is because, as shown in [4], the binding energy per particle of larger nuggets increases monotonically as a function of N, asymptoting to the saturation value at very large N. Up to angularmomentum-dependent effects, it is thus energetically favorable to form successively larger nuggets, with the threshold for dissociation of nuggets growing ever larger than the temperature. 2 By arguing through analogy with the SM, in the next subsection we discuss the circumstances under which a nugget bottleneck could appear.
A feature of large nuggets discussed in [4] is that the total binding energy per dark (baryon) number can be a substantial fraction of m X . Thus if very large nuggets are synthesized, because all of this binding energy is released in the synthesis process, one might worry that substantial heating of the dark bath relative to the SM bath could occur. However, under reasonable assumptions, we will see that the binding energy release results in at most an Oð1Þ change in the dark sector temperature. Assuming that the dark bath consists of only the scalar mediators, ϕ, and the dark matter, we may place an upper bound on the change in temperature: where BE Nsat =N is the saturation binding energy per dark number. 3 For , the temperature change will be 2 This is in contrast to large nuclei in the SM, where binding energy per nucleon grows with mass number A up until a global maximum near A ∼ 60; the drop is due to Coulomb (repulsion) energy and symmetry energy associated with neutron-proton imbalance. 3 If all of the dark matter were bound in nuggets of effectively infinite dark number and there were no additional components of the dark bath, then the inequality would be saturated. negligible. The analogous statement for nucleosynthesis is T BBN ∼ 100, which, of course, is easily satisfied. Here since ϕ is massive and because the binding energy per dark number can be an order-one fraction of m X , we can does not necessarily hold. However, since the asymmetric abundance n X has frozen out well before synthesis begins, 5 heat release increases the ϕ number density but not the X number density, which serves to limit the overall temperature change. Taking the equilibrium number density X in the parameter range of interest. The temperature change can approach order 1 when ϕ is nonrelativistic at the onset of synthesis (when BE 2 ∼ m ϕ , corresponding to m ϕ =m X at the upper range of what we consider) or when ϕ is very highly relativistic m ϕ ⋘ m X at the onset of synthesis. In the relativistic ϕ limit, wherem X is the average mass per dark number and η is the baryon-to-photon ratio. When m ϕ ⋘ m X , thenm X ¼ ] assuming most DM is bound into saturated nuggets at the end of synthesis. Though n X =n ϕ can be larger than naively expected, as long as ffiffiffiffiffiffiffiffiffiffiffiffi ffi Putting this all together in the m ϕ ⋘ m X limit we then have As discussed above, BE Nsat =N in the m ϕ ⋘ m X limit. So e.g. for α ϕ ∼ 0.1, one could approach an order-one temperature change when ffiffiffiffiffiffiffiffiffiffiffiffi ffi m ϕ m X p ≲ 10 −4 GeV or an order of magnitude when ffiffiffiffiffiffiffiffiffiffiffiffi ffi m ϕ m X p ≲ 10 −8 GeV. It is important to note that the heat per dark number, ∼BE Nsat =N, is not released all at once, but rather in increments over the entire synthesis process up to large nuggets. Since on average the dissociation temperature rises as a function of nugget size, the heat release should not delay synthesis.
C. Bottlenecks at low dark nugget number?
By arguing through analogy with the SM, here we discuss the circumstances under which a dark nugget bottleneck could appear. First we briefly review the nature of the SM BBN bottlenecks. BBN proceeds when the temperature drops below the deuteron binding energy (order MeV) so that deuteron dissociation is suppressed. Once a substantial fraction of deuterons exists, almost all free neutrons are incorporated into 4 He [with intermediate steps of e.g. DðD; nÞ 3 He, 3 HeðD; HÞ 4 He]. A few factors prevent any significant amount of synthesis to higher-A isotopes. First, the Coulomb barrier is substantial at T BBN and becomes rapidly more substantial as A and Z grow. Second, building up from 4 He to larger nuclei through Coulomb-barrier-free neutron capture is stifled because 5 He, with its unpaired neutron, is unstable to decay to 4 He þ n. Additionally, a 4 He þ 4 He fusion process is endothermic, including 8 Be þ γ, such that 8 Be is unstable. Taken together, these facts imply the presence of bottlenecks not only at A ¼ 2, but also at A ¼ 5 and A ¼ 8. This guarantees a small abundance of nuclei larger than 4 He exiting BBN.
We can consider the existence of bottlenecks in nugget synthesis by analogy with the SM. By construction there is no obstruction due to a Coulomb barrier. However, fluctuations in binding energy as a function of nugget size, leading to unstable states analogous to A ¼ 5, 8 for nuclear matter, could lead to substantial bottlenecks. The absence of stable A ¼ 5, 8 nuclei can be understood qualitatively within the shell model of the nucleus. Within the shell model, constituents are treated as noninteracting, but each is under the influence of an emergent potential. Constituents fill up available states from low to high energy, obeying Pauli's exclusion principle. The existence of many approximately degenerate energy eigenstates results in some states having significantly larger binding energy per constituent than their neighbors with slightly larger nugget number. The constituent numbers N i of these exceptional states are referred to as magic numbers. For spherically symmetric bound states, the first two magic numbers are always 2 and 8 (corresponding to doubly magic 4 He and 16 O nuclei), corresponding to the filling of l ¼ 0 and 1 states. For larger N, the locations for the magic numbers depend on the potential and spin-orbit coupling, and may scale like N i ∼ i 2 to i 3 (for Coulomb and quadratic potential).
States with size just over or equal to a sum of magic numbers are prone to instability or metastability. Isospin 4 However, recall that α ϕ ≳ 0.01. See Fig. 1. 5 Since T syn ≲ BE 2 =15 ∼ α 2 ϕ m X =60, X is nonrelativistic at the onset of synthesis when the coupling is perturbative. During synthesis, the average mass per dark number may decrease as more X's are bound into (more and more tightly bound) nuggets, but any individual DM nugget will always have mass greater than or equal to m X . slightly complicates the story for nuclei, but consider states closely related to doubly magic 4 He. The doubled state, 8 Be, decays to 2 4 He. States with one extra neutron or proton ( 5 He or 5 Li) decay to 4 He plus a free nucleon. For the A ¼ 5 states, pairing also plays a role; nuclei with unpaired nucleons of given isospin are less strongly bound. This is because (pseudoscalar) pions mediate an attractive interaction between opposite-spin nucleons, effectively reducing the energy of the configuration. Since we have only a scalar mediator, it is unclear whether an analog of a pairing effect might exist. 6 In any case, in the context of small nuggets, if the binding of 4 X is smaller than twice the binding energy of 2 X, the decay process 4 X → 2 X þ 2 X will become possible. Similarly, perhaps the rest energy of 3 X is greater than that of 2 X þ X since the third nucleon in 3 X sits alone, unpaired, in an l ¼ 1 shell. For larger N, the degeneracy of energy eigenstates is typically broken more strongly and the spacing of energies becomes smaller so that magic numbers become rare and unimportant.
Since the system under consideration contains only one type of constituent (isospin zero) instead of two as in the SM, as just discussed, the analogue to unstable A ¼ 5 and A ¼ 8 states for nuclei is unstable N ¼ 3 and N ¼ 4 nuggets. If both 3 X and 4 X are unstable and have small lifetimes in comparison to the cross section for fusions ðn 3ð4Þ vσð 3ð4Þ X þ 2 X → 5ð6Þ X þ 0 XÞÞ −1 , then there is a bottleneck to fusion of larger nuggets. In the absence of a Coulomb barrier, if only one of 3 X and 4 X were unstable, it would be easier to synthesize past this artificial bottleneck (e.g. 3 X þ 2 X → 5 X þ ϕ or 2 X þ 2 X → 4 X þ ϕ). 7 Since the shell model is fundamentally phenomenological and not predictive in terms of a UV completion, we cannot make definitive statements about the presence or absence of these bottlenecks. We will thus frame our discussion of nugget synthesis in terms of the presence or absence of bottlenecks.

IV. NUGGET SYNTHESIS
We have just reviewed the conditions for the formation of an N ¼ 2 nugget ( 2 X). Once a population of 2 X exists, there are a number of routes to creating larger nuggets. Here we discuss synthesis of larger-N nuggets given either the presence or absence of low-N bottlenecks. First we review fusion processes and reasonable models for their respective cross sections. Next we set up the Boltzmann equations that govern the evolution of the mass function and introduce a dimensionless interaction time-a dimensionless function of dark number density, typical cross section, typical speed, and Hubble rate-that sets the nugget size scale at the end of synthesis. Then we provide analytic estimates of typical nugget size after synthesis in both the presence and absence of bottlenecks. Finally we present semianalytic solutions to the Boltzmann equations for the mass function after synthesis and relate them to our analytic estimates. Our semianalytic discussion of the solution to the Boltzmann equations is similar to that of [7], though with two significant differences. First, utilizing the results of our companion paper, we are able to compute the synthesized nugget mass in terms of the parameters of the model, m X , m ϕ , and α ϕ . Second, we employ the compound nucleus model, and account for the possibility of reactions resulting in final states that include an extra nugget emission, which can lead to a broader mass function exiting synthesis.

A. Dominant processes and cross sections
The relevant reactions, ignoring any contribution from N-body interactions for N > 2, are of the form, where 0 X denotes the mediator ϕ, and X number conservation dictates that i þ j ¼ k þ l þ m þ Á Á Á Here we describe and motivate using the compound nucleus (CN) model of nuclear physics for modeling cross sections for these processes. Within the CN model, one treats such reactions in two stages: first, an excited CN of size i þ j is formed, and then the compound state disintegrates [12] (see also e.g. [13]). If any excess energy within the CN is quickly shared amongst all constituents, then it is reasonable to treat the two stages as independent processes. In this case, the probability for disintegrating into any given final state depends only on the energy, angular momentum, and parity of the CN, with no specific dependence on how the state is formed. We can expect this so-called Bohr assumption to apply in most nugget reactions involving at least one saturated nugget in the early Universe. This is because the characteristic time for energy to be transferred across a nugget of size R is order R=v s where v s ¼ dp=dϵ is the speed of sound within a nugget, which we find to be a substantial fraction of the speed of light for saturated nuggets. On the other hand, the time scale over which the nuggets interact is order R ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 1 − b=R p =v, where b is the impact parameter and v the relative speed of the initial state nuggets. For large nuggets, v is typically much smaller than the speed of light; thus we expect the excess energy to be efficiently distributed within the compound system over the time of a typical interaction as long as b ≲ R. There is a potential exception coming from interactions with b ≃ R, where grazing interactions could lead to two-to-two dark number exchange processes.
Given the Bohr assumption, the cross section for reactions Eq. (13) takes the form 6 By virtue of paired constituents being in close proximity, since real scalars will mediate a strong attractive interaction independent of spin, one can imagine that a similar pairing effect might arise. 7 Though note that rates for these processes could be low if where σ ij is the cross section for forming the compound state of size i þ j and Γ α;iþj is the rate for the compound state to disintegrate into a final state, α. We will now address models for the formation cross section σ ij and the disintegration rates, Γ α;iþj , in turn. For collisions of two large, saturated nuggets, we adopt a geometric cross section, with r 0 given by Eq. (3). For collisions between a large saturated nugget and a much smaller nugget (i ≫ j) we adopt a quantum-corrected geometric cross section similar to models used for neutron capture, q is the small nugget momentum and is the effective momentum inside the large nugget, and m eff ¼ jμ is the effective mass of the small nugget inside i X. The 1=p correction to the radius accounts for the effective de Broglie wavelength of the small nugget, which leads to an enhancement for small nugget capture. For larger nuggets sizes R i , 1=p is typically subdominant and thus neglected. The transmission factor T accounts for quantum reflection effects due to a sudden change in the small nugget's effective mass upon entering the large saturated nugget; note that in the limit that the smaller nugget is also saturated, p ¼ p 0 and thus T ¼ 1.
Time reversal invariance or reciprocity relates disintegration rates to formation cross sections for the reverse process. Consider the decay of a CN, N X Ã , to a less excited state N−k X Ã , while emitting a much smaller nugget k X (or a mediator for k ¼ 0). It can be shown (see Refs. [13,14]), when the density of states into which the CN can decay is large, that the partial width of the CN to decay into products N−k X Ã þ k X is described by a thermal distribution weighted by the energy release [using the liquid drop model in Eq. (2)] in the decay, where ΔE Ã ≈ ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi M 2 k þ p 2 q − kμ is the energy release. TðE Ã Þ is the temperature of the nucleons in the excited CN, − M k is the kinetic energy of the emitted particle, g k the degrees of freedom of the emitted particle, and σ k is the cross section for the inverse process N−k X Ã þ k X → N X Ã . For very low temperatures satisfying T ≪ M k − kμ, it is immediately apparent that emission of the lightest possible particles is heavily favored if g k σ k is not radically different for different k. In our case, we will see that the temperatures are sufficiently small that we expect ϕ emissions to dominate.
At low temperatures, the nugget is described as a Fermi gas with a modified fermion mass m Ã . Then, to leading order in the excitation energy E Ã , the temperature is is the Fermi energy given by Eq. (3). The binding energy release in the process, fN X þ ð1−fÞ X → N X, is given by Thus, assuming negligible kinetic energy in comparison to binding energy release, the maximal excitation energy of a CN is E Ã ∼ 0.26ϵ surf N 2=3 . This leads to maximum temperatures of order where in the final expression we used the definition N sat ≡ ðr 0 m ϕ Þ −3 and Eq. (3). In Ref. [4], we found that ϵ surf is order 3m X to 10m X for sample parameters α ϕ ¼ 0.1, 0.01 and m ϕ =m X ¼ 10 −2 ; 10 −3 ; 10 −4 , and it mildly increases with decreasing m ϕ =m X and decreasing α ϕ . At saturation, the typical temperature of the formed CN will be much smaller than m X but possibly comparable to or larger than m ϕ . In the parameter ranges of interest (where saturation applies), scalar emissions generally dominate by many orders of magnitude over small nugget emission, with stronger domination as N grows and m ϕ =m X falls.
Despite domination of scalar emission over small nugget emission at each step in the cascade, one might worry that the large number of steps in a CN cascade could lead to a sizable total small nugget emission rate. An upper bound on the number of X emissions can be estimated as N X ≲ N ϕ Γ X ðE Ã 0 Þ=Γ ϕ ðE Ã 0 Þ with N ϕ the number of ϕ emissions given an initial excitation energy E Ã 0 . Using the further upper-bound estimate N ϕ ≲ E Ã 0 =m ϕ , we find that N X is substantially below 1 for the relevant parameter ranges discussed in Sec. III A.

B. Boltzmann equations
We have just seen in Sec. IVA that coagulation-fusion through emission of mediators only-will typically dominate when the Bohr assumption of quick thermalization applies. At larger impact parameter, a thermal model may not apply and nugget fragmentation could occur. We will suppose that-if they are relevant at all-the dominant form of nonthermal fusion reactions are of the form i X þ j X → iþj−k X þ k>0 X þ ϕ's, which we will refer to as two-to-two reactions.
In either case, the nugget distribution will follow a set of Boltzmann equations. Defining y k ≡ nk X =n X with n k the number density of a nugget of size k and n X ¼ P k kn k , we have and for notational simplicity, we denote coagulation as a two-to-two process with a final state particle 0 X. All fusion processes will generally include many additional ϕ emissions, which do not impact the nugget number distribution. At low temperature, only exothermic processes contribute, and the summations in Eq. (18) are restricted to final states that are more asymmetric than the initial state. Total dark number conservation implies that P k ky k ¼ 1, so that in the large N continuum, ky k becomes a probability distribution. There are two sources of time dependence for the Boltzmann equation: the density n X , which dilutes away as ∼1=a 3 , and the thermal-averaged cross sections, which may include a nontrivial transmission factor T . It is possible to factor out the time dependence by defining a dimensionless interaction time γ, such that (cf. w in [15]) where σ ∘ ¼ πr 2 0 , and hvT i ∘ is the common factor obtained from a thermal average of the velocity-dependent part of the cross sections. For fusion of nuggets both deep in the saturation regime, T ¼ 1, and we simply take hvT i ∘ ≃ ffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi T X =m X p . Note that the binding energy per particle can be very large and thusm X can be significantly smaller than m X . For 2 X-large fusion, relevant in bottleneck scenarios, we instead take T ∼ 4v, and hvT i ∘ ≃ 4T X =m X . In terms of γ, the Boltzmann equation becomes where the factor σ ∘ hvT i ∘ absorbs the time dependence of hσvi, such that the temporal evolution of the nugget distribution is entirely captured by γðtÞ.
In general, the nugget bath temperature T X can deviate from the temperature of the photon bath T γ , if the dark sector is not in kinetic equilibrium with the SM. Depending on whether the nugget bath is relativistic (T X ∝ a −1 ) or not (T X ∝ a −2 ), the form of γðtÞ varies. In the deep saturation limit during radiation domination, one has and for the case of 2 X-large nugget fusion, one has where the subscript "syn" denotes the value of the quantity at the beginning of synthesis. For BE 2 not radically greater than m ϕ , ϕ is nonrelativistic at T syn already, and without an additional dark relativistic component, the dark bath will be nonrelativistic.
As discussed later, the typical nugget size exiting synthesis will be of order γ 6=5 max with γ determined through Eq. (21) in the absence of a bottleneck, and the typical size of a subdominant population of very large nuggets will be order γ 3 max with γ determined through Eq. (22) when a strong bottleneck is present.

C. Large nugget formation in the absence of a bottleneck
Here we derive an analytic understanding of the typical size of a nugget exiting synthesis. We start by defining an average nugget size. Motivated by the fact that P k ky k ¼ 1, which indicates ky k acts like a probability distribution in the continuum limit, we define the average size bȳ k ≡ P k¼1 k 2 y k . 8 For large k, m k ∝ km X , andk is essentially an energy density-weighted average. The evolution ofk follows Here the summation only includes processes where the final states have a more asymmetric nugget distribution, 8 An alternative definitionk ¼ P k¼1 ky k P k¼1 y k ¼ ð P k y k Þ −1 corresponding to a number density-weighted average is not appropriate when two-to-two processes contribute significantly toward fusion of large nuggets. This is because P k y k does not change under two-to-two processes, but only under fusion through mediator emissions. such that i þ j − k > j ≥ i > k, implying thatk is monotonically increasing.
To proceed further, we take the saturation limit, which will be a good approximation as long ask ≳ N sat . Consider the most dominant contributions in Eq. (23), which comes from fusions of large nuggets with sizes aroundk, and emissions (mostly) of mediator particles ϕ. This process roughly doubles the size, and leads to an Oð1Þ change in logk. When the time scale for this process to occur exceeds γ max , synthesis freezes out. Therefore, we approximate Eq. (23) by replacing the cross section (summed over k) in Eq. (23) with the totalk −k interaction rate, and substitute ðk − jÞðk − iÞ ∼k 2 and P i;j y i y j ∼ 1=k 2 . The freeze-out condition then becomes In the saturation limit, the cross section fork −k fusion scales like σ¯k¯k ∼ σ ∘k 2=3 , with a velocity dependence v¯k¯k ∼ v ∘k −1=2 . (The transmission factor is T ¼ 1 in this limit.) Then, setting hσvð¯kX þ¯kX → ∼2kÞ Xi ¼ σ ∘ v ∘k 1=6 , the typical nugget size at the end of synthesis is with γ max given by Eq. (21) when a syn =a → 0, and the typical nugget mass isM fo ∼k fomX . Explicitly, where fiducial values correspond roughly to the parameters α ϕ ¼ 0.1, m X ¼ 10 GeV and m ϕ ¼ 10 MeV. In Fig. 2 we show contours of constantk fo andM fo =GeV assuming saturation values for r 0 andm X as described in [4] and summarized in Sec. II B, as a function of m X and m ϕ , with two choices of fixed α ϕ . We have taken T syn γ ¼ T syn X , with T syn X estimated as described in Sec. III B, and we have included the contribution of ϕ to g Ã when relevant. 9 The m X range shown is that which satisfies the synthesis condition Eq. (6) and the condition T syn X > 10T eq . The shaded region at smaller m X =m ϕ does not satisfy the synthesis condition in Eq. (5). We chose the lower cutoff for m ϕ according to ffiffiffiffiffiffiffiffiffiffiffiffi ffi m ϕ m X p ≳ 10 −8 GeV. For smaller values of ffiffiffiffiffiffiffiffiffiffiffiffi ffi m ϕ m X p , as discussed in Sec. III B, if the SM and DM sectors are kinetically decoupled and the thermal DM sector does not contain additional light degrees of freedom, the dark sector temperature could increase by an order of magnitude or more during synthesis, which would increase the estimate by a similar factor.
In Fig. 2, we can see thatk fo andM fo depend much more strongly on m ϕ and α ϕ than on m X . Using BE 2 ¼ α 2 ϕ m X =4 and the estimates Eq. (3), modulo the very weak FIG. 2. Contours of typical nugget number exiting big bang darkleosynthesis,k fo (dashed red line) and typical nugget massM fo (solid purple) for α ϕ ¼ 0.03 (left) and α ϕ ¼ 0.3 (right). The temperature of the dark sector is assumed to be roughly T X ≈ T γ . The blue shaded region corresponds to BE 2 < m ϕ , where 2 X synthesis will not be efficient [Eq. (5)]. The upper m X cutoff corresponds to the requirement that the 2 X fusion rate be smaller than Hubble as in Eq. (6), and the lower m X cutoff corresponds to T synth ≲ 10T eq . (See Fig. 1.) The various kinks in the contours are results of the change in g Ã as the synthesis temperature passes through the QCD phase transition and neutrino decoupling. 9 As discussed in Sec. III B, depending on details of the cosmology, we expect T syn γ ¼ T syn X to within a factor of a few, depending on the decoupling time of the two sectors. dependence of g Ã and BE 2 =T syn X on model parameters, we find that the typical nugget numberk fo and massM fo scale as This scaling is reflected in the plot. For example, it is worth noting thatk fo grows slightly more rapidly with decreasing m ϕ than doesM fo . This is entirely due to the increase of binding energy per particle with decreasing m ϕ . Recall that our estimate is valid only whenk fo > N sat since we have used the saturation cross section and binding energy in our freeze-out estimate. In the strong binding limit, Thus sincek fo scales more strongly with m −1 ϕ , for fixed α ϕ and m X , the approximation becomes better as m ϕ decreases. Numerically, we find that for α ϕ ≲ 0.1, the self-consistency condition for applying the saturation limit k fo > N sat is always satisfied whenever the synthesis conditions are met. In Fig. 2, where the region with k fo < N sat for α ϕ ¼ 0.3 is shaded red, we see that even for larger α ϕ the estimate is invalid only in a small region of parameter space.

D. Nugget distribution function in the absence of a bottleneck
So far we have focused on order of magnitude analytic estimates. Here we support and complement our estimates with a full semianalytic analysis of the Boltzmann equations in Eq. (20), complemented with numerical calculations. As before, we only consider coagulation and two-to-two processes. However, much of the analysis carries over when higher order processes are included, as long as they stay homogeneous with the same weight as we will discuss. It is useful to rewrite the infinite set of differential equations in Eq. (20) by introducing a kernel, Kði; j; kÞ,  [17]. Here we consider the saturation limit and utilize the CN-like picture for two-to-two processes, such that the kernel scales simply as Kði; j; kÞ ¼ where Γ k is proportional to the partial width of a compound state i þ j transitioning into a final state k þ ði þ j − kÞ, the squared factor characterizes the scaling of the geometric cross section, and the square root factor characterizes the relative speed. A similar kernel was considered in [7], but with Γ k ¼ δ iþj;k , corresponding to the case of coagulation [17]. 10 There are generally no closed form solutions even for a simplified choice of fusion kernel [16], and given that k and the total interaction time γ max are typically very large, numerical calculations quickly become intractable. Fortunately, as similarly discussed in [7], when the kernel is homogeneous Kði; j; kÞ ¼ λ −α Kðλi; λj; λkÞ, scaling solutions exist, which allow for extrapolation of numerical results for small γ to large γ ∼ γ max . In the saturation limit, the total inelastic cross section is already homogeneous with degree 1=6. Then a scaling solution is valid as long as Γ k in Eq. (31) is also homogeneous, which we will assume for now. It is important to note that the kernel, Eq. (31), applies only when large saturated nuggets dominate the dark matter density; for fusion involving at least one unsaturated nugget, the form of the kernel will change, and for the most general fusion interactions, one cannot even define a time-independent kernel. We will check that our results are self-consistent in that the vast majority of nuggets exiting synthesis are above the saturation size, where the kernel we use does apply.
We will now derive the scaling solution, following [16] in spirit. In the largek limit, nugget indices may be treated as continuous variables. Consider the ansatz y k ðγÞ ¼ s 2 fðksÞ, where sðγÞ is some function of γ that is to be determined. Substituting the ansatz into Eq. (30) and replacing summation with integration in the continuum limit, we have 10 Specifically, Ref. [7] used Kði;jÞ ¼ ði −1=2 þ j −1=2 Þði 2=3 þ j 2=3 Þ to match with existing mathematical literature on the coagulation equations. However, as we explicitly show, it is not necessary to choose such a kernel for the scaling solution to apply. Therefore, we use the more physical kernel that properly characterizes the total geometric cross section and relative velocity. where we have used the homogeneity property of K to change the integration variable. We have not explicitly included the integration bounds, which do not affect the derivation as long as they are linear functions of the integration variables. One sees that for s ¼ cγ 1=α , the γ dependence drops out entirely, and one is left with an integro-differential equation for fðxÞ, given by For nugget fusion, we will consider K with homogeneous weight α ¼ −5=6, and hence s ∼ γ −6=5 . Equation (33) in general is still very difficult to solve even numerically. However, it is known that the nugget distribution generally approaches the scaling solution very quickly independent of the initial condition (see [16]). Therefore, it is possible to numerically integrate Eq. (30) by truncating the differential equation at finite nugget number, and then extract the scaling function fðxÞ by testing that the solutions y k ðγÞ have converged to a scaling limit. It is illuminating to revisit our earlier discussion of typical nugget sizek in Sec. IV C in light of the scaling solution. In the scaling limit,k ¼ s −1 R dx x 2 fðxÞ. Given that there is freedom to choose fðxÞ by rescaling s, we may setk ¼ s −1 . Then the scaling limit simply becomes k 2 y k ðγÞ → k 2 k 2 f k k : ð34Þ In the scaling limit, k 2 y k maintains the same shape, centered onk, with a scaling behaviork ∼ γ 6=5 , verifying our earlier estimate. For our numerical study, we consider three separate homogeneous Γ k for the kernel Kði; j; kÞ. The three branching ratio forms we consider correspond to (i) Coagulation: Where Γ ij→k ∼ δ iþj;k (ii) Uniform: Γ ij→kl ∼ θðjk − lj − ji − jjÞ.
(iii) Energy scaling: Γ ij→kl ∼ Q 2 ij→kl σ kþl , where the heat release Q ij→kl is proportional to the change in total surface area for a given reaction. Specifically, which roughly captures the increase in phase space as the reaction becomes more exothermic. Coagulation dominates when (a) the compound nucleus model applies, (b) the difference between the binding energy per particle of very large and very small nuggets is greater than the mediator mass, and (c) the excited nucleus temperature is small [see Eq. (17) and surrounding discussion]. In Sec. IVA we argued that the compound nucleus model should apply except possibly when the impact parameter is comparable to the nugget radii-i.e. for grazing collisions. While points (b) and (c) are valid for the parameter ranges we consider, it is useful to have a sense for what might change in models where they do not hold. Therefore we consider potential corrections to the coagulation-only picture by analyzing two other reasonable choices for the branching ratios where our scaling solution still applies.
Starting with the initial condition y i ð0Þ ¼ δ i0 , Fig. 3 shows k 2 y k ðγÞ as a function of nugget number k and a specific nugget branching fraction as a function of interaction time γ. In the continuum limit, k 2 y k ðγÞ is interpreted as the differential nugget probability distribution, dpðkÞ=d log k, or the FIG. 3. Nugget distribution function at different interaction time γ (right) and specific number fraction as a function of γ. Here, the different line style depicts different assumptions for the branching ratio, with the solid/dashed/dotted lines corresponding to the energyscaling/uniform/coagulation branching ratio assumption as described in the main text. fraction of the nugget number shared by nuggets centered on k within a unit bin in log k. We see that the nugget distributions very quickly become dominated by large k even at small interaction time γ ¼ 100. Differences in the branching ratios do not significantly alter the behavior, although a flatter branching ratio tends to enhance the nugget distribution at small k. Figure 4 shows the extracted mass function at different interaction time γ ∈ f10; 20; 100g, in dotted/dashed/solid line for the different branching ratio assumptions. One already sees a convergence to a universal function. At smaller nugget number, the mass functions are significantly broader for branching ratios that include two-to-two processes. However, the difference between the uniform and energy scaling branching ratios are small.

E. Large nugget formation in the presence of bottleneck: Nugget capture
In this section we consider nugget synthesis in the presence of a bottleneck. As discussed in Sec. IV E, we expect the nugget distributions to be impacted significantly only when both 3 X and 4 X are unstable. In this case, the two-to-two fusion processes are halted until larger nuggets ( 5 X or 6 X) are formed through other processes. For instance, 6 X can be produced via rare three-2 X fusions, where a shortlived state 4 X may exist, permitting the reaction 4 X þ 2 X → 6 X þ 0 X to proceed before the 4 X decays. This is analogous to 3 4 He → 12 C through the 8 Be resonance. In this case, the 6 X now functions as a nucleation site, capturing nearby 2 X.
In the limit that the nucleation sites are sparse, the process N X þ 2 X → Nþ2 X with ϕ emissions dominates. Considering only this coagulation process, for large k, the Boltzmann equation in Eq. (20) becomes Working in the saturation limit, and taking the large k continuum limit, the cross sections will scale as hσvð k X þ 2 X → kþ2 XÞi ≃ σ ∘ hvT i ∘ k where fiducial values correspond roughly to the benchmark parameters α ϕ ¼ 0.1, m X ¼ 10 GeV, m ϕ ¼ 10 MeV. The wave equation description breaks down when 2 X starts to be depleted, which happens when the fractional DM number contained in the nucleation sites becomes Oð1Þ. If we assume that the probability of a single 2 X squeezing through the bottleneck at the beginning of synthesis is 2p, then the number density of the nucleation sites will be pn X . The nucleation sites will evolve linearly until the fractional DM number in nucleation sites, pðξγÞ 3 , is roughly 1=2, at the transition interaction time At this point, the average nucleation size will bē k ∼ 1=ð2pÞ. Beyond this point, the 2 X distribution is expected to rapidly become depleted, and the subdominant large-large fusion will become significant. If γ 3 trans < γ 6=5 max , the discussion in Sec. IV C will then apply once again, with a scaling lawk ∼ γ 5=6 . Figure 5 shows the nugget sizes and masses for the nucleation sites exiting synthesis. We have assumed that p is small enough that the 2 X population always dominates. Compared to Fig. 2, the final nugget number and masses are significantly larger due to γ 3 max scaling. This is expected as the fusion rate is controlled by the 2 X densities which remain relatively large. One may be concerned that the local 2 X density within a Hubble volume of a nugget may be depleted before γ max is reached, which would render our calculation invalid. Such a requirement will impose a lower bound on p ≳ H 3 =n X ∼m X GeV T 3 ηm 3 Pl , which is negligibly small when the temperature is of order GeV or smaller.

F. Nugget distribution and bottlenecks
Here we present a numerical investigation of the effects of bottlenecks on the nugget distribution. In Sec. IV E we argued that a significant bottleneck might occur due to 3 X and 4 X both being unstable, while isolated unstable nuggets do not cause significant changes in the nugget distributions.
In order to simulate nugget evolution in this situation, we follow the same setup in Sec. IV D, and use the cross sections in Eq. (31). Though Eq. (31) is modeled on largelarge fusion, the velocity dependence approaches a constant for small-large fusions and describes the bottleneck scenario well. The large-large fusion cross section is not fully captured by the kernel, as its temperature dependence is different form the small-large case. However, contributions from large-large interactions are subdominant whenever , which is relevant for the benchmark cases we have. The Boltzmann equation is truncated at k ≤ 1000, such that all branchings into higher k nuggets are fixed to zero. We also utilize the energy scaling branching ratios, though different branching assumptions do not change the quantitative behavior significantly. The initial conditions are chosen such that 2y 2 ¼ 1 − p and 6y 6 ¼ p. To simulate the effects of bottlenecks, we fix y k ¼ 0 and Γ kl ¼ 0 if k X or l X is unstable. The left panel of Fig. 6 shows the mass distributions for different choices of p. We see that in general, the nugget distributions are separated into two distinct populations: the small 2 X population and large nugget nucleation sites. The distributions for the nucleation sites are strongly peaked, and move forward rapidly as γ increases. The total fraction of DM in the nucleation sites increases as well. Variations in p simply modify the total nucleation site density. At large γ, one expects the 2 X density to appreciably decrease eventually, although a FIG. 5. Contours of the nugget number, k fo (dashed red line), and the typical mass of the nuggets, M fo (solid purple) for the large nugget populations exiting darkleosynthesis when a significant bottleneck is present. The coupling is fixed at α ϕ ¼ 0.03 (left) and α ϕ ¼ 0.3 (right). Similar to Fig. 2, the blue shaded region corresponds to inefficient 2 X fusion due to small binding energy, and the upper (lower) cutoff for m X corresponds to an 2 X fusion rate always smaller than Hubble and T synth ≲ 10T eq respectively. numerical confirmation is impractical, as it requires including exponentially more terms in the Boltzmann equations.
The right of Fig. 6 also shows a comparison of different nugget functions when including extra bottlenecks at k ¼ 9 for the dashed line (and also 10 for the dotted line), while fixing p ¼ 10 −5 . With only one extra bottleneck at k ¼ 9, the mass functions quickly move beyond the bottleneck and become indistinguishable from the ones without the k ¼ 9 bottleneck. While for the case with two bottlenecks at k ¼ 9, 10, no larger nuggets are produced beyond k ¼ 8. This confirms our expectations, where isolated bottlenecks do not change the nugget distribution while multiple consecutive bottlenecks can bring fusion processes to a halt.

V. CONCLUSIONS
We have studied the early Universe synthesis of manyparticle bound states of ADM, or nuggets. We focused on a minimal model with spin-1=2 asymmetric dark matter interacting through a scalar mediator; the scalar is solely responsible for binding the nuggets and determines the dynamics of fusion in both the initial and final stages of synthesis. We unified the treatment of synthesis with a quantitative calculation, in our companion paper [4], of the properties of the large bound states, utilizing effective field theory tools developed in nuclear physics. Within our model, the typical nugget size exiting early Universe synthesis can be many orders of magnitude larger than has been estimated in previous treatments of many-particle bound states of dark matter. We derived a nugget mass function, describing the dark matter energy distribution in nuggets of different sizes, from the Boltzmann equation. We argued that within the minimal model fusion occurs primarily through radiation of mediators, and showed that in this case the mass function is peaked within a couple of generations in size about the typical size, confirming the findings in [7]. We found that inclusion of the analog of baryon-number-transfer nuclear fusion reactions in addition to radiative fusion reactions generates a broader nugget distribution.
With a quantitatively derived nugget mass function exiting early Universe synthesis, we are now in a position to study the late Universe cosmology of ADM nuggets. This is the subject of our next investigation.