Implications of an Improved Neutron-Antineutron Oscillation Search for Baryogenesis: A Minimal Effective Theory Analysis

Future neutron-antineutron ($n$-$\bar n$) oscillation experiments, such as at the European Spallation Source and the Deep Underground Neutrino Experiment, aim to find first evidence of baryon number violation. We investigate implications of an improved $n$-$\bar n$ oscillation search for baryogenesis via interactions of $n$-$\bar n$ mediators, parameterized by an effective field theory (EFT). We find that even in a minimal EFT setup, there is overlap between the parameter space probed by $n$-$\bar n$ oscillation and the region that can realize the observed baryon asymmetry of the universe. The mass scales of exotic new particles are in the TeV-PeV regime, inaccessible at the LHC or its envisioned upgrades. Given the innumerable high energy theories that can match to, or resemble, the minimal EFT that we discuss, future $n$-$\bar n$ oscillation experiments could probe many viable theories of baryogenesis beyond the reach of other experiments.

Introduction -The search for physics beyond the Standard Model (BSM) requires efforts at both high energy and intensity frontiers. In this regard, a particularly powerful probe is offered by rare processes that violate (approximate) symmetries of the Standard Model (SM), such as baryon and lepton numbers (B and L), which can be inaccessible to high energy colliders but within reach of low-energy experiments. A well-known example is proton decay, whose non-observation leads to strong constraints on ∆B = ∆L = ±1 new physics even at the scale of Grand Unified Theories (GUTs), ∼ 10 16 GeV [1,2].
Baryon and lepton number violation are intricately tied to one of the outstanding puzzles in fundamental physics, the origin of the baryon asymmetry in the universe. If baryogenesis occurs at temperatures above the weak scale, B − L violation is required to avoid washout by electroweak sphalerons. In this regard, constraints from proton decay (which conserves B − L) are not applicable. Here we consider instead B-violating, Lconserving new physics at an intermediate (sub-GUT) scale, so that baryogenesis may proceed both above and below weak scale temperatures. From the low energy point of view, effects of heavy new particles are encoded in higher dimensional operators in an effective field theory (EFT), where B-violating, L-conserving interactions can appear first at the dimension-nine level [3], in the form of |∆B| = 2, ∆L = 0 operators. In this case, neutron-antineutron (n-n) oscillation (see [4] for a recent review) is well placed to search for B violating phenomena and shed light on baryogenesis. 1 1 Other |∆B| = 2, ∆L = 0 processes include dinucleon decays: nn → π 0 π 0 , pp → π + π + , pn → π + π 0 probe the same operators as n-n oscillation, but with a lower sensitivity at present [5], while pp → K + K + [6] can be relevant for B-violating new physics with suppressed couplings to first-generation quarks [7,8].
Current measurements constrain the free neutron oscillation time to be τ nn > ∼ 10 8 s [9,10]. Upcoming experiments, in particular at the European Spallation Source (ESS) and also potentially the Deep Underground Neutrino Experiment (DUNE), are poised to improve the reach up to 10 9-10 s [11][12][13][14][15]. As we will see in detail below, such numbers translate into new physics scales of roughly (τ nn Λ 6 QCD ) 1/5 ∼ O(10 5-6 GeV), well above the energies directly accessible at existing or proposed colliders. Discussions of the physics implications of a potential n-n oscillation discovery, in particular for baryogenesis, are therefore both important and timely.
We now elucidate the connection between τ nn and the new physics scale in the EFT context. The lowest dimension effective operators contributing to n-n oscillation at tree level are dimension-nine operators of the form O nn ∼ (uudddd). The classification of these operators dates back to the 1980s [54][55][56][57][58] and was refined recently in [59], which established an alternative basis more convenient for renormalization group (RG) running. A concise review of the full set of tree-level n-n oscillation operators is provided in the Appendix. In what follows, we focus on one of these operators for illustration, Here u, d are SM up and down quark fields, respectively, and u c , d c are their charge conjugates. i ( ) , j ( ) , k ( ) are color indices, and "h.c." denotes hermitian conjugate. The operator suppression scale Λ (1) nn is generally a weighted (geometric) average of new particle masses, modulo appropriate powers of couplings and loop factors.
If the operator is generated by integrating out new particles at a high scale M , computing τ nn requires RG evolving the EFT down to a low scale µ 0 (usually chosen to be 2 GeV), where it can be matched onto lattice QCD. The leading contribution to RG rescaling reads [58,59]  Here α (n f ) s is the effective strong coupling with n f light quark flavors, whose value is obtained with the RunDec package [60]. Corrections from two-loop running as well as one-loop matching onto lattice QCD operators were recently computed [59] and are small, and will be neglected in our calculations. No additional operators relevant for n-n oscillation are generated from RG evolution.
The n →n transition rate is determined by the matrix element of the low-energy effective Hamiltonian between the neutron and antineutron states. Thus, once n|O nn (µ 0 )|n are known, we can relate τ nn = n|H eff |n −1 to the six-quark operator coefficients. Recent progress in lattice calculations [61,62] has greatly improved the accuracy and precision on n|O nn (µ 0 )|n compared to previous bag model calculations [56,57] often used in the literature. Using the results in [62], and assuming the operator in Eq. (1) gives the dominant contribution to n-n oscillation, we can translate the Super-K limit into Λ (1) nn > ∼ 4 × 10 5 GeV (for a representative RG rescaling factor of 0.7). An improvement on τ nn up to 10 9 (10 10 , 10 11 ) s will correspond to probing Λ (1) nn ∼ 5 (8, 13) × 10 5 GeV. These numbers are representative of the whole set of n-n oscillation operators, and do not vary significantly with the starting point of RG evolution M (see Appendix for details).
A minimal EFT for n-n oscillation and baryogenesis -One of the simplest possibilities for generating the operator in Eq. (1) at tree level is with a Majorana fermion X of mass M that couples to the SM via a dimensionsix operator of the form 1 Λ 2 Xudd, which originates at an even higher scale Λ M via some UV completion that we remain agnostic about. A familiar scenario that realizes this EFT setup is supersymmetry (SUSY) with R-parity violation (RPV), where the bino plays the role of X and the dimension-six operator is obtained by integrating out squarks at a heavier scale. However, this simple EFT with a single BSM state does not allow for sufficient baryogenesis due to unitarity relations: in the absence of B-conserving decay channels, X decay cannot generate a baryon asymmetry at leading order in the B-violating coupling, a result known as the Nanopoulos-Weinberg theorem [63] (see [64] for a recent discussion); meanwhile, 2 → 2 processes uX →dd andūX → dd are forced to have the same rate and thus do not violate CP .
A minimal extension that can accommodate both n−n oscillation and the observed baryon asymmetry involves two Majorana fermions X 1 , X 2 (with M X1 < M X2 ), each having a B violating interaction 1 Λ 2 Xudd. In addition, a B conserving coupling between the two is necessary to evade constraints from unitarity relations. In the context of RPV SUSY, this corresponds to the presence of a wino or gluino in addition to the bino, which is known to allow for sufficient baryogenesis [64][65][66].
Guided by minimality, we assume X 1,2 are both SM singlets, and consider just one of the many possible B conserving operators in addition to the two B violating ones. Our minimal EFT thus consists of the following dimension-six operators that couple X 1,2 to the SM: 2 2 Our minimal EFT bears similarities with the models studied in [67,68]. However, these papers focused on baryogenesis using operators of the form (d c P R d)(ū c P R X), which, upon Fierz transformations, are equivalent to generation-antisymmetric components of the (ū c P R d)(d c P R X) operators in Eq. (3), and thus do not mediate n-n oscillation at tree level.
Both X 1 and X 2 mediate n-n oscillation -integrating them out at tree level gives This setup contains all the necessary ingredients for baryogenesis [69]: the Lagrangian in Eq. (3) violates B and P , while nonzero phases of η X1 , η X2 , and η c can lead to CP violation; departure from equilibrium can occur in multiple ways, as we discuss below. Although a clear simplification, we expect the minimal set of operators in Eq. (3) to capture the generic qualitative features possible in a two n-n mediators setup, which can be realized in more complicated and realistic frameworks.
We calculate the baryon asymmetry by numerically solving a set of coupled Boltzmann equations to track the abundances of X 1,2 and B − L (B) above (below) T = 140 GeV (we assume sphalerons are active when T > 140 GeV, resulting in Y B = 28 79 Y B−L ). Our aim is to find regions of parameter space that can achieve the observed Y B = 8.6 × 10 −11 [70,71], with suitable choice of CP phases. Technical details of this calculation can be found in the Appendix.
If all three operator coefficients have similar sizes, Λ X1 ∼ Λ X2 ∼ Λ c , it is difficult to obtain the observed baryon asymmetry in the region of parameter space probed by n-n oscillation. For M X1,2 > ∼ 10 4 GeV, the Λ's that can be probed are sufficiently low for X 1,2 to remain close to equilibrium until their abundances become negligible, while efficient washout suppresses B(−L) generation. For lower masses and higher Λ's, on the other hand, X 2 may freeze out with a significant abundance, and decay out of equilibrium at later times when washout has become inefficient, so that both limitations from the Late decay

Early decay
FIG. 1. Sketches of the evolution of the heavier n-n mediator abundance YX 2 , washout rate Γwo and baryon asymmetry YB in the two scenarios considered in this letter (arbitrary normalization). In the late decay scenario, the n-n mediator is long-lived and decays out of equilibrium to generate a baryon asymmetry. In the early decay scenario, departure from equilibrium (thin dotted curve) is small, but suppressed washout enables efficient baryogenesis. See text for details.
higher mass regime are overcome. However, its CP violating branching fraction Achieving the desired baryon asymmetry in the ESS/DUNE reach region therefore requires hierarchical Λ's; such scenarios can arise if new particles in the UV theory that mediate the corresponding operators have hierarchical masses and/or couplings, or if the EFT operators are generated at different loop orders. We find compatible regions of parameter space in two distinct scenarios, one with late decays of X 2 and the other with earlier decays. These are schematically illustrated in Fig. 1, and discussed in turn below (a detailed analysis with benchmark numerical solutions is presented in the Appendix).
Late decay scenario -For Λ X2 ∼ Λ c Λ X1 , nn oscillation is dominated by X 1 exchange and probes the M X1 -Λ X1 parameter space (see Fig. 2). This hierarchy leads to weaker interactions for X 2 compared to the degenerate case, causing it to freeze out with a higher abundance Y fo X2 . Also, X 2 becomes long-lived and decays after washout processes have become ineffective, thereby creating substantial baryon asymmetry (see Fig. 1). In this case, its CP -violating branching fraction scales as and does not decouple as Λ X2 and Λ c are both increased, enabling Y B ∼ Y fo X2 CP to reach the observed value. Numerically, we find that this baryogenesis scenario is viable with Λ X2 , Λ c > ∼ 20 Λ X1 in the parameter space probed by n-n oscillation. In Fig. 2, we show regions in the M X1 -Λ X1 plane that can accommodate the observed Parameter space of the minimal EFT probed by n-n oscillation for the late decay scenario, assuming MX 2 = 4 MX 1 . For ΛX 2 = Λc = 50 ΛX 1 , the green shaded region can accommodate YB = 8.6 × 10 −11 . For ΛX 2 = Λc = 25 ΛX 1 (100 ΛX 1 ), viable region is between dashed red (dot-dashed blue) lines. The gray shaded region marks ΛX 1 < MX 2 , where EFT validity requires greater than O(1) coupling.
baryon asymmetry for various choices of Λ X2 /Λ X1 = Λ c /Λ X1 . In each case, the lower boundary of the viable region is effectively determined by the requirement that X 2 freezes out with sufficient abundance. As we move upward from this lower boundary, increasing all three Λ's while keeping their ratios fixed, at some point we enter a regime where X 2 decouples from the SM bath while relativistic, and Y fo X2 saturates at Y eq s , so that further increasing the Λ's only reduces CP and hence the final Y B . Furthermore, for sufficiently high Λ X2 and Λ c , X 2 dominates the energy density of the universe before it decays (this does not happen for X 1 in the parameter space we consider), so that its decay injects significant entropy into the plasma, diluting the baryon asymmetry. Both of these effects -saturation and dilution -determine the upper boundary of the viable region.
Early decay scenario -For the opposite hierarchy Λ X1 Λ X2 , n-n oscillation is dominated by X 2 exchange and probes the M X2 -Λ X2 parameter space (see Fig. 3). In this case, X 2 is short-lived, and its abundance closely follows the equilibrium curve. However, small departures from equilibrium, always present in an expanding universe because interaction rates are finite, can be sufficient for baryogenesis if washout can be suppressed. The rates for washout processes involving X 1 and X 2 are proportional to n 1 Λ −4 X1 and n 2 Λ −4 X2 , respectively, where n 1,2 are the number densities of Parameter space of the minimal EFT probed by n-n oscillation for the early decay scenario, assuming MX 2 = 4 MX 1 . Points represent solutions with YB = 8.6 × 10 −11 found in a scan over ΛX 2 < ΛX 1 < 100 ΛX 2 , MX 2 < Λc < ΛX 2 . For all these points, ΛX 1 ∼ 10 ΛX 2 is needed to suppress washout. The gray shaded region marks ΛX 2 < MX 2 , where EFT validity requires greater than O(1) coupling.
would be efficient until T ∼ M X1 , i.e. until n 1 starts to fall exponentially. In contrast, by increasing Λ X1 , we enter a regime where washout is dominated by X 2 processes at high temperatures and becomes inefficient as soon as the temperature falls below M X2 (washout due to udd ↔ūdd, whose rate ∼ T 11 /M 2 Λ 8 falls steeply with T , is also irrelevant at this point), resulting in a short period of baryon asymmetry generation from X 2 decays (see Fig. 1). Note that increasing Λ X1 with respect to Λ X2 also helps to increase departures from equilibrium compared to the degenerate case. Fig. 3 shows points in the M X2 -Λ X2 plane that can realize the observed Y B through this early decay process, based on a numerical scan over the region Λ X2 < Λ X1 < 100 Λ X2 , M X2 < Λ c < Λ X2 . For the majority of these points, Λ X1 is within a factor of two from 10 Λ X2 , while Λ c < ∼ 3 M X2 . The results can be understood from the competing effects of baryon asymmetry generation and washout, where the rate of baryon asymmetry generation Γ ∆B =0 is calculated from CP -violating X 2 decays. First of all, a lower ratio Λ c /M X2 is always preferable (within the range of EFT validity), while the ratio Λ X2 /Λ X1 has an optimal value of ∼ 1/10 as a result of balancing between faster baryon asymmetry generation at higher tempera-tures (which favors higher Λ X2 /Λ X1 ) and later transition to X 1 -dominated washout (which favors lower Λ X2 /Λ X1 ). The requirement of sufficient departure from equilibrium precludes arbitrarily low Λ c and leads to a minimum M X2 for this scenario to work, which we see from Fig. 3 is a few ×10 4 GeV. Finally, the overall size of Λ X1,2 is essentially determined by the requirement that Y B freezes out around the time Γ ∆B =0 X2 /Γ wo reaches its maximum, and is higher for higher M X2 .
Complementary probes -In the region of parameter space that is allowed by existing n-n oscillation searches, within reach of the ESS/DUNE, and realizes the observed baryon asymmetry, we find M X1,2 > ∼ 10 3 (10 4 ) GeV in the late (early) decay scenario. Given that X 1,2 are SM singlets that only couple to the SM via higher dimensional operators, it is unlikely that they can be detected at the LHC or its envisioned upgrades. Likewise, there are no strong flavor physics constraints on our minimal EFT with just the operators in Eq. (3). We note, however, that this outlook can change in a more complicated model that preserves the general features of baryogenesis of our minimal EFT if at least one of X 1,2 carries SM charges or couples to other fermion species. For example, colored particles at the TeV scale, such as the gluino in RPV SUSY, could be within LHC reach. Likewise, extending the exotic fermion couplings to other quark flavors can introduce potential constraints from flavor violation considerations such as K 0 −K 0 mixing [67]. Nevertheless, our minimal EFT study illustrates that n-n oscillation might be uniquely placed to probe realistic baryogenesis scenarios that are inaccessible via other searches.
Conclusions -Establishing baryon number violation (or the absence thereof up to a certain scale) will have far-reaching implications on our understanding of fundamental particle interactions, in particular on the mechanism that generates the observed baryon asymmetry in our universe. Motivated by the unprecedented sensitivity to n-n oscillation that can be achieved at future facilities, the ESS and DUNE in particular, which offers new opportunities to probe |∆B| = 2, ∆L = 0 interactions, we studied implications of a potential discovery for baryogenesis scenarios involving n-n mediators. We took a bottom-up EFT approach with a minimal set of fourfermion operators coupling the n-n mediators to the SM, which, despite being simplistic, sets a useful template that more sophisticated theories can build upon. We identified two viable baryogenesis scenarios -one involving late out-of-equilibrium decays of a heavy Majorana fermion, and another involving earlier decays assisted by a suppressed washout rate -that can be realized in the parameter space to be probed by future n-n oscillation searches, with no corresponding collider or flavor signals. These results highlight the capability of n-n oscillation experiments to probe an important BSM phenomenon, that of baryogenesis, beyond the scope of other searches.
We thank J. Barrow, G. Brooijmans and C. Csáki for helpful comments and discussions. C.G. was supported by the European Commission through the Marie Curie Career Integration Grant 631962, and by the Helmholtz Association. B.S. was partially supported by the NSF CAREER grant PHY-1654502 and thanks the CERN and DESY theory groups, where part of this work was conducted, for hospitality. The work of J.D.W. was supported by the DoE grant de-sc0007859 and the Humboldt Research Award. The work of Z.Z. was supported by the DoE grant de-sc0007859, the Rackham Dissertation Fellowship, and the Summer Leinweber Research Award. J.D.W. and Z.Z. also thank the DESY theory group for hospitality. This work was performed in part at the Aspen Center for Physics, which was supported by National Science Foundation grant PHY-1066293.
Here we briefly review the effective operator analysis of n-n oscillation. Since multiple operators may be present in addition to the representative operator we considered in the letter, to gain intuition about the new physics scale being probed, let us define (A.1) As we will see explicitly below, Λ nn defined here roughly coincides with suppression scales of dimension-nine operators mediating n-n oscillation. This is because the nuclear matrix elements n|O nn |n ∼ O(Λ 6 QCD ). Taking  where the number in Eq. (A.2) shows the current best limit from Super-K. There are 12 independent operators that contribute to n-n oscillation at tree level. Using the basis of [59], we write andŌ i is obtained by exchanging P L ↔ P R in O i . Note that since QCD conserves parity, O i andŌ i have identical nuclear matrix elements and anomalous dimensions (neglecting weak interactions). Our labeling of O 1,2,3 is in accordance with [59], while our O 4,5,6 are proportional to their Q 5,6,7 , respectively (their Q 4 , which we have skipped here, has zero nuclear matrix element). The operator basis of Eq. (A.5) is particularly convenient because different operators do not mix as they are evolved from some high scale(s) µ (i) down to µ 0 = 2 GeV, where lattice calculations of nuclear matrix elements are reported. We have Here we have chosen µ = 10 5 GeV as a reference scale to compute the numbers, and introduced r i ≡ α s (µ i )/α s (10 5 GeV),r i ≡ α s (μ i )/α s (10 5 GeV) to account for effects due to different choices (when O i andŌ i are renormalized at µ i andμ i , respectively, rather than at 10 5 GeV).
In the special case that the RHS of Eq. (A.7) is dominated by a single term, say the one proportional to c i ≡ Λ we can establish a correspondence between τ nn (equivalently Λ nn ) and Λ nn . This is shown in Fig. 4. As mentioned above, all Λ

Boltzmann equations
The Boltzmann equations to be solved for our minimal EFT are dn a dt + 3Hn a = C a (a = 1, 2, 3) , (A.8) where n 1,2 are the number densities of X 1,2 , and n 3 represents n B−L (n B ) for T > 140 GeV (T < 140 GeV) when electroweak sphalerons are assumed to be active (inactive). We define where f eq a is the equilibrium distribution at zero chemical potential for species a. Assuming a common temperature is maintained for all species, we have f a = e µa/T f eq a ≡ r a f eq a ≡ (1 + ∆ a ) f eq a , (A.10) for the actual distribution of species a, with ∆ a characterizing the amount of departure from equilibrium. The collision terms can then be written in terms of the W 's and r's, where W udd→ūdd , W ūdd→udd are computed from the corresponding matrix elements with contributions from onshell X 1,2 exchange subtracted. We have grouped together terms that are identical as dictated by CP T invariance, W i→f = Wf →ī (where bar denotes CP conjugate state).
To further simplify, we note that several processes conserve CP up to one-loop level, and as a result For the CP -violating processes, on the other hand, we define their CP -symmetric and CP -asymmetric components, As a consequence of CP T invariance and unitarity, f W i→f = f Wī →f , which implies Using these relations and noting that r d rd = 1 (because µ d + µd = 0), the collision terms can be rewritten as As |∆ u,d | 1 in all cases, we have only kept terms up to linear order in ∆ u,d . In addition, we approximate ∆ u,d = e µ u,d /T − 1 µ u,d /T . We have dropped the W (0) udd→ūdd term, which is higher order in 1/Λ. Now the collision terms are written in terms of ∆ 1,2,u,d , while the LHS of the Boltzmann equations contain n 1,2,B(−L) . To relate the two sets of quantities, we note that, assuming Maxwell-Boltzmann distributions for X 1,2 , ∆ a = n a n eq Meanwhile, the chemical potentials µ u,d are related to n B(−L) (see e.g. [72]): for T > 140 GeV, as follows from equilibration of Yukawa interactions and SU (3) and SU (2) sphalerons, and conservation of hypercharge; for T < 140 GeV, dT . This is the final form of the Boltzmann equations that we use in our numerical solutions. In order to determine viable parameter space regions for baryogenesis, we set arg(η * X1 η X2 η c ) = π/2 to maximize CP violation, and look for solutions with the final Y B ≥ 8.6 × 10 −11 ; for such parameter choices, the exact amount of observed baryon asymmetry can then be achieved with some suitable choice of arg(η * X1 η X2 η c ) ≤ π/2. As mentioned in the letter, if X 2 is sufficiently long-lived, its decay may dump significant entropy into the plasma, diluting the baryon asymmetry. We account for this effect by dividing the final Y B from solving the Boltzmann equations by a dilution factor d s = 1.83 h Here Γ X2 is the total decay width of X 2 , and x d and Y d X2 are the values of M X2 /T and Y X2 at the time of X 2 decay, determined by Γ X2 = H.

Interaction rates
We now provide analytical expressions for the interaction rates W that appear in the collision terms. For a 2 → 2 process ab → cd, W ab→cd = n eq a n eq b σv ab→cd = where S i , S f are symmetry factors for the initial and final states (e.g. S i = 2 if a and b are identical particles) and in the center of mass frame. The sum is over initial and final state spins and colors, while " " means averaging over cos θ, with θ being the scattering angle in the center of mass frame. We take the upper limit of integration to ∞ for simplicity since the integrand is exponentially suppressed for center-of-mass energies above the EFT cutoff Λ for temperatures where the EFT is valid (T Λ). Computing the scattering amplitudes at tree level, we find where φ c = arg η c . We have seen above that all the CP violation in 2 → 2 processes can be encoded in W uX1→uX2 . Computing also one-loop diagrams for this process, we find We have explicitly checked that CP violation in uX a →dd satisfy expectations from unitarity relations Eq. (A. 16). For decay processes, where Γ is the rest frame decay rate, summed over final state spins and colors, and averaged over the initial state spin. At tree level, we have where ρ = M X1 /M X2 . The CP -violating decay rate at one-loop level reads This function is maximized at ρ = 0.265 1/4.

Benchmark solutions
To have a more detailed understanding of the baryogenesis scenarios discussed in the letter, let us examine a few benchmark solutions to the Boltzmann equations. We choose M X2 = 4M X1 = 2 × 10 5 GeV, which can accommodate solutions in both the late and the early decay scenarios, and consider the following three benchmarks: • Degenerate: Λ X1 = Λ X2 = Λ c = 1.5 × 10 6 GeV.
All three benchmarks induce n-n oscillation at a level that is consistent with current constraints and may be within reach of future searches.
We plot the evolution of various quantities from solving the Boltzmann equations in Fig. 5. The upper-left panel shows the amount of departure from equilibrium for X 1 (dashed) and X 2 (solid), quantified by ∆ a = (Y a − Y eq a )/Y eq a , while the solid curves in the upper-right panel show the baryon asymmetry Y B .
We first note that, with the exception of ∆ 2 in the late decay benchmark, departures from equilibrium are always very small due to efficient depletion of X 1,2 number densities by rapid decays once they become nonrelativistic. As a rough estimate, assuming radiation domination, we have Γ 1→3 /H ∼ 10 −5 M 5 when T < M , where 10 −5 is the size of the phase space factor. For M ∼ 10 5 GeV and Λ ∼ 10 6 GeV, we have Γ 1→3 /H > ∼ 10 5 and thus efficient decays that keep ∆ a 1. On the other hand, ∆ 2 in the late decay benchmark evades this pattern with much higher values for Λ 2,u ∼ 10 8 GeV, which result in Γ 1→3 /H ∼ 10 −3 (M X2 /T ) 2 , and thus later decay, for X 2 . In this case, ∆ 2 starts to grow exponentially once the most efficient 2 → 2 process dX 2 →ūd freezes out, which happens when Γ 2→2 /H ∼ n eq d σv /H ∼ g eff Setting them equal, we have This power law dependence is clearly seen from the upper-left panel of Fig. 5. Also note that ∆ a is larger for higher Λ, as it is harder to catch up with Hubble expansion when interactions are weaker. When nondegenerate Λ's are involved in the X 1,2 number changing processes, the lowest of them tends to determine the total interaction rate, and thus ∆ a . For example, at high temperatures, ∆ 1,2 are lower in the early decay benchmark compared to the degenerate case because of a lower Λ c . They exceed the degenerate curves later when coannihilation becomes Boltzmann suppressed; from here on, ∆ 1 tends to grow faster due to a higher Λ X1 , while the lower Λ c maintains ∆ 1 ∆ 2 via u X 1 ↔ u X 2 , X 2 ↔ X 1ū u processes. Next, to understand the trend of the Y B curves in the upper-right panel of Fig. 5, it is useful to note that the Boltzmann equation for Y B has the following form,