Dark Matter, Shared Asymmetries, and Galactic Gamma Ray Signals

We introduce a novel dark matter scenario where the visible sector and the dark sector share a common asymmetry. The two sectors are connected through an unstable mediator with baryon number one, allowing the standard model baryon asymmetry to be shared with dark matter via semi-annihilation. The present-day abundance of dark matter is then set by thermal freeze-out of this semi-annihilation process, yielding an asymmetric version of the WIMP miracle as well as promising signals for indirect detection experiments. As a proof of concept, we find a viable region of parameter space consistent with the observed Fermi excess of GeV gamma rays from the galactic center.


Introduction
The existence of dark matter (DM) has been confirmed through many independent studies, and we now know that this non-baryonic matter comprises around a fifth of the energy density of the universe [1][2][3]. If DM is composed of a single new particle, then the Weakly Interacting Massive Particle (WIMP) paradigm is particularly attractive, since the presentday abundance of DM can be set by thermal freeze-out in an expanding universe (see [4] for a review). Given that the standard model (SM) features a vast and subtle structure, though, there may be comparably rich dynamics in the dark sector. Perhaps the dark sector features some of the same properties as visible SM baryons, such as accidental symmetries, asymmetric abundances, and instability of some of its components.
In this paper, we propose a new DM scenario where asymmetries in the visible and dark sectors are closely connected. Due to light unstable states carrying baryon number, asymmetry sharing is efficient down to temperatures below the DM mass, such that the ultimate (asymmetric) abundance of DM is set by thermal freeze-out of a semi-annihilation process. This can be regarded as a hybrid framework between asymmetric and WIMP scenarios, whereby the processes responsible for asymmetry sharing in the early universe can potentially produce signals today in indirect detection experiments.
As a prototypical example, consider a stable DM particle A which carries baryon number 1/2. This DM particle couples to an unstable particle B with baryon number 1 and to a  Figure 1. Two approaches to sharing an asymmetry between the DM and SM sectors: (a) through a higher-dimension operator that becomes inefficient as the universe cools, or (b) though an on-shell intermediate particle B which remains in quasi-equilibrium until A freezes out. In this paper, we assume that B chemically decouples from the SM before A freezes out (indicated by the dashing), in which case freeze-out of semi-annihilation AA → Bφ controls the final DM abundance.
light mediator φ with baryon number 0. The semi-annihilation [5][6][7][8][9][10] process allows the SM baryon asymmetry to be shared in the early universe, and it also gives rise to promising indirect detection signals in the galactic halo today. Note that this process does not involve A, since this is an asymmetric DM scenario. As we will see, there is a viable region of parameter space that not only yields the desired DM and baryon abundances but also yields a GeV gamma ray signal consistent with the galactic center (GC) excess seen in Fermi [11][12][13][14][15][16][17][18][19][20][21][22]. The presence of asymmetry sharing is familiar from asymmetric DM models (see [23,24] for reviews), but there is a key difference highlighted in figure 1. 1 In scenarios shown in figure 1(a), some high-dimension operator or heavy state links the visible and dark sectors [30] (see also [31][32][33][34][35][36][37][38][39][40]). As the temperature of the universe drops, this interaction becomes irrelevant, locking in independent asymmetries for baryons and DM. In our scenario shown in figure 1(b), the interactions that link the visible and dark sectors involve light mediators in quasi-equilibrium, and therefore stay relevant to energies below the DM mass. As long as m B < 2m A , eq. (1.1) is always kinematically allowed, impacting both thermal freeze-out in the early universe and indirect detection today.
There are three relevant temperatures in this scenario. The decoupling temperature T D is when the mediator B chemically decouples from SM baryons, setting the asymmetries in the dark sector. The effective symmetry imbalance temperature T I is when the abundance of A (and B) departs significantly from the abundance of A (and B). The freeze-out temperature T F is when the semi-annihilation process AA → Bφ freezes out, setting the final A abundance. We will focus on the hierarchy which yields an asymmetric version of the WIMP miracle. As we will see, achieving

JCAP02(2016)052
With this mass range it is natural (though not required) to consider B as a bound state of two A particles. For example, the semi-annihilation process in eq. (1.1) could correspond to the formation of dark bound states [41], dark nuclei [42,43], or dark atoms [44], though the key difference compared to previous work is that the resulting B particle is unstable. 2 Unlike standard asymmetric DM models, eq. (1.1) leads to indirect detection signals without requiring a residual symmetric component. The resulting indirect detection spectrum depends on the decays of B (and to a lesser extent, the decays of φ). As benchmarks, we consider two possibilities involving SM quarks, 3 B → udd, B → cbb, (1.4) and both of these decay modes can yield a good fit to the Fermi GC excess, albeit with different preferred values of m A and m B . As we will see, the desired phenomenology requires m A and m B to have comparable masses of O(10-100 GeV). 4 As a proof of concept, we highlight a benchmark scenario that satisfies the following three criteria, summarized in figure 2.
• Correct DM abundance: We require Ω DM h 2 ∈ [0.115, 0.124], within the 95% Planck confidence limit of the DM abundance [3]. 5 • Correct baryon asymmetry: We require η b = (n b − n b )/s ∈ [8.1, 9.5] × 10 −11 , within the 95% confidence limit of the baryon asymmetry of the universe [46]. 6 • Plausible fit to the GC excess: We follow the strategy of ref. [22] to find the best fit masses for the GC excess gamma ray spectrum.
Simultaneously satisfying these criteria turns out to be possible, since much of the intuition derived from symmetric WIMPs holds also in this asymmetric scenario. We check that this benchmark is consistent with CMB heating bounds from Planck [3], and we evaluate possible direct detection and antiproton flux constraints. We also investigate searches for displaced jets at the Large Hadron Collider (LHC), which could be a distinctive signature for the proposed scenario in optimistic regions of parameter space. The rest of this paper is organized as follows. We introduce our prototype model in section 2 and study its cosmological history and asymmetry sharing in section 3, leaving details to the appendices. In section 4, we show how the AA → Bφ semi-annihilation process followed by B decay can match the observed GC excess gamma ray spectrum. We discuss constraints and potential collider signals in section 5 and conclude in section 6.
2 To achieve a more predictive framework, one could consider B to be a bound state of A particles mediated by φ exchange. In this case, mB would be determined by mA, m φ , and the A-φ coupling, and there would be tight relationship between the various interaction cross sections. In order to explore the full parameter space, we will not assume that B is a bound state resulting from φ interactions. Indeed, in dark nuclei scenarios [42,43], B is bound because of (residual) confining interactions and φ is a spectator to the binding process. 3 These kinds of decay modes also appeared in the asymmetric DM scenario of ref. [45]. In that model, B itself was DM and had late-time decays. Here, A is DM and B decays are prompt on cosmological scales. 4 The astute reader may notice that the above processes do not conserve fermion number. This is easily reconciled with an extended model presented in appendix A, yielding equivalent phenomenology. Similarly, φ can correspond to multiple unstable states, specifically two in appendix A. 5 For simplicity, we assume a Gaussian distribution to extrapolate the 2 σ Planck range from the 1 σ limit (ΩDMh 2 ∈ [0.1175, 0.1219]) in ref. [3]. 6 Here, s is the entropy density. This differs from the more familiar notation for η b , which is defined with respect to the photon density. Today s = 7.05nγ [46,47].

Abundances & GC Best Fit
Example region of parameter space where we obtain the desired DM abundance [3] (red), the desired baryon asymmetry [46] (blue), and a good fit to the GC excess from AA → Bφ semiannihilation followed by B → udd (purple) or B → cbb (orange) decays. Here, Ω DM h 2 and η b are plotted with respect to the true cross section σ AA→Bφ v , while the GC fit regions are shown with respect to the effective cross section σ AA→Bφ v eff scaled by the DM abundance (see eq. (4.5)). The hashed region indicates the CMB heating bounds from eq. (5.1), evaluated with respect to σv eff . Not shown are the AMS-02 antiproton bounds from eq. (5.5), which do constrain this parameter space but have large uncertainties. The white star indicates the benchmark parameters in eq. (2.7).

Prototype of Dark Matter with a shared asymmetry
The analysis in this paper is based on a simplified dark sector consisting of three particles: a stable DM species A with baryon number 1/2, an unstable state B with baryon number 1, and a light mediator φ which couples to the SM through, e.g., the hypercharge portal [48][49][50][51][52], Higgs portal [53,54], or axion portal [55,56]. The reason the A particle is exactly stable is that there are no states in the SM with baryon number 1/2. An explicit Lagrangian with these properties is presented in appendix A, where the state A is replaced with a fermion/boson system. For simplicity, we take this ψ A /φ A system be mass degenerate in our discussion, such that the DM dynamics can be captured by an effective single species A. It is also possible to split ψ A from φ A to achieve a more varied phenomenology.
There is considerable freedom in choosing the masses of A and B, though we focus on scales relevant for describing the GC excess: The mass of the φ is assumed to be small compared to the other scales in the theory, m φ O(10 MeV-1 GeV), consistent with φ being a light dark photon that mixes with the SM photon or a light scalar that mixes with the Higgs sector. We take the couplings of φ to SM states to be large enough for φ to stay in thermal equilibrium in the early universe, but small enough to avoid direct detection bounds on A (see section 5). The key interactions in this scenario are shown in figures 3 and 4. The processes in figure 3 control the freeze-out abundance of A, and for simplicity, we assume these processes have velocity-independent cross sections. The semi-annihilation process in figure 3(a) has various crossed and conjugate channels of relevance: In order to have an asymmetric DM model with a suppressed symmetric component, we assume that the annihilation processes in figures 3(b) and 3(c), are highly efficient. When needed for numerical studies, we take the amplitudes |M AA→φφ | = |M BB→φφ | = 1, though smaller couplings lead to similar results. 7 In the extreme asymmetric limit with negligible A and B abundances, the only (thermally-averaged) cross section of relevance is σ AA→Bφ v . 8 The processes in figure 4 control the sharing of the SM baryon asymmetry with the dark sector. We denote the total asymmetry by η tot , which will later be divided into the baryon asymmetry η b and DM asymmetry η A . Assuming B is a fermion with decay modes given by eq. (1.4), the B decays in figure 4(a) are mediated by the operators

4)
7 For such large annihilation amplitudes, one might wonder whether φ exchange could result in additional bound state formation between A and B particles (see, e.g., refs. [57][58][59]). Such bound states could be avoided either by reducing the φ coupling strength (perhaps resulting in a small residual symmetric DM component) or by ensuring that the would-be Bohr radius of the bound state is larger than the φ Compton wavelength (see appendix A). Alternatively, the dominant annihilation process could be unrelated to the φ particle, as expected if A and B were part of a larger strongly-coupled dark sector with both stable and unstable components. 8 When φ is sufficiently light, the (semi-)annihiliation process can be Sommerfeld enhanced and therefore velocity dependent. This adds non-trivial temperature dependence to the evolution of the comoving A and B abundances and can change the viable region of parameter space. We neglect the Sommerfeld effect in this work, though one generically expects it to boost the present day semi-annihiliation rate compared to the rate during DM freezeout.  where Λ O(10-10 3 TeV) is some new physics scale. 9 The B decay width is given parametrically by (2.5) This short lifetime ensures that B decays prior to the start of big bang nucleosynthesis (BBN), T 10 MeV, corresponding to t 10 −2 sec [46,61]. The same operators also contribute to the 2 → 2 scattering processes involving B shown in figure 4(b), which is the mechanism responsible for the initial sharing of the asymmetry between the SM and the DM in the early universe (see figure 5 below).
As we will see in the next section, there are three temperatures of interest: the B chemical decoupling temperature T D , the A/A (and B/B) symmetry imbalance temperature T I , and the A freeze-out temperature T F . While T D depends sensitively on the scale Λ, as long as T D > max{T I , T F } then the details of B decoupling are irrelevant for the evolution of the Boltzmann system. Furthermore, if T I T F , then the details of the A annihilation are relatively unimportant. Assuming the hierarchy T D > T I T F , the dominant phenomenology of this scenario can be determined by four parameters: (2.6) Below, we refer to the following benchmark parameters, which yields the desired behavior specified in the introduction with the B → cbb decay channel: Note that this benchmark cross section value is comparable to canonical WIMP scenarios. The proximity of m A and m B is because the GC excess fit in section 4 prefers a higher Lorentz boost factor for B. 10 As shown in section 3.4, we get reasonable DM phenomenology whenever m A m B < 2m A . 9 If B is a scalar, then the decay must include an additional fermion, perhaps a SM lepton or a supersymmetric particle (see, e.g., ref. [60]). If B is a fermionic bound state of a fermion/boson ψA/φA pair, then one should regard eq. (2.4) as effectively a dimension 7 operator of the form ψAφAu c d c d c , such that the effective cutoff scale Λ would be lower than our benchmark value. 10 In addition, as shown in figure 7(b), taking the ratio mB/mA to be in the vicinity of 1 gives fine control over the DM abundance for fixed semi-annihilation cross section, making it easier to simultaneously fit the GC excess and achieve the desired DM abundance. One might be surprised that η tot in this benchmark is comparable to the present-day baryon asymmetry η b 9 × 10 −11 , such that DM only takes around 1/20 of the asymmetry. This is unlike usual asymmetric DM scenarios where the total asymmetry η tot is roughly twice the baryon asymmetry η b [23,24]. The main reason for the mismatch is that m A is around ten times heavier than the standard 5 GeV benchmark for asymmetric DM (i.e. the mass that yields the right DM abundance when η DM η b ). Also note that the concrete model in appendix A has both a boson and fermion species of A, and A has a baryon number 1/2 (see eq. (3.1) below), further affecting the expected mass relation.

Initial conditions
We start from a primordial asymmetry generated through an unspecified mechanism (see [23,24] for a review of options). The total asymmetry η tot is then conserved through the thermal history of the universe by the processes in figures 3 and 4. We use the notation Y i = n i /s, where n i is the number density of species i and s is the entropy density. We can express the total asymmetry as where the factor of 1/2 accounts for the baryon number of A and η b is the baryon asymmetry (baryons minus anti-baryons). For ease of discussion, we refer to DM as a single particle A. Because B is a fermion, however, our explicit model in appendix A requires A to be a degenerate fermion-boson system, a fact which is reflected in the equations below. At high temperatures, the SM and DM asymmetries are related since the operators in eq. (2.4) (or their ultraviolet completions) ensure that the two sectors are in chemical equilibrium. A full analysis of the chemical potentials is presented in appendix B. In equilibrium at a temperature T , the asymmetries are: where for the explicit model in appendix A, The factor of 45/29 is familiar from ordinary asymmetric DM scenarios, see e.g. [62]). The function f is defined in eq. (B.10); it has the asymptotic behavior f f (x) → 1/6 for fermions and f b (x) → 1/3 for bosons as x → 0 (early times) and f (x) → 0 as x → ∞ (late times).

Chemical decoupling of B
The equilibrium conditions in eqs. 10 15 Reaction Rates Figure 5. Comparison of the scattering rate of Bq → qq and the decay rate of B → qqq with the Hubble expansion rate. By eq. (3.5), the decoupling temperature T D occurs when the scattering rate falls below the expansion rate, which in turn sets the initial conditions for the subsequent Boltzmann evolution of A particles. Chemical equilibrium between B particles and the SM is restored when the B decay rate becomes relevant, though in our benchmark studies, this happens after A particles have already frozen out.
the B particles chemically decouple from the SM. To estimate the decoupling temperature T D , we compare the rate of B scattering/decay to the Hubble expansion where the thermally-averaged rate for Bq → qq scattering is calculated in appendix C, σv Bq→qq ∼ T 2 /Λ 4 in eq. (C.14), the number density of quarks n eq q0 ∼ T 3 is given in eq. (C.9), the chemical potential of quarks µ q is defined in eq. (B.18), and the B decay width Γ B is given in eq. (2.5). The Hubble parameter is where g * is the number of relativistic degrees of freedom at the temperature T [47]. Because B interacts with the SM via the contact operators in eq. (2.4), the scattering process are more relevant at early times (high temperatures) while the decay process is more relevant at late times (low temperatures). 11 In figure 5, we compare the B scattering/decay rates with the Hubble expansion rate for various choices of Λ, taking the benchmark dark sector parameters in eq. (2.7). We plot temperatures in terms of x = m A /T . For Λ 50 TeV, decoupling never happens since B and the SM remain in chemical equilibrium throughout the early history of the universe. This tends to erode the DM asymmetry η A , since for T m A , f (m A /T ) → 0 in eq. (3.2). Eventually, the semi-annihilation process AA → Bφ freezes out, which stops the depletion of η A . While this could lead to a potentially interesting phenomenology, in this paper we focus on larger values of Λ where B decouples prior to A freeze out. Taking the benchmark JCAP02(2016)052   Figure 6. Evolution of the comoving abundances of the {A, B} system, using (a) the benchmark parameters from eq. (2.7) and (b) a more bound-state-like spectrum with m B /m A = 1.9. The dashed curves indicate the equilibrium distributions, and the dotted curves indicate the anti-particle distributions. The arrow shows the approximate location of the semi-annihilation freezeout y F . In (b), it is coincidental that the decay of B starts around the same time as y F . in eq. (2.7) as an example, Λ = 300 TeV yields T D 55 GeV, which corresponds to x 1.1 in figure 5.
When the dark sector decouples from the SM at T D , the asymmetries can be estimated by evaluating eqs.
, and the asymmetries at the time of decoupling are: For the benchmark in eq. (2.7) with T D 55 GeV, the actual values are which is not so different from the T D → ∞ limit.

Thermal freeze-out of A
Once B chemically decouples from the SM, the dark sector follows standard Boltzmann evolution. The full Boltzmann system is described in appendix D, which includes the impact of having initial asymmetries for A and B. There are two relevant dimensionless time variables, where the numerator of y was chosen to match the approximate kinetic energy available in the AA → Bφ process. Note that x > y whenever m B < 2m A . In figure 6(a), we show the numerical evolution of the various Y X for the benchmark scenario in eq. (2.7). As x increases, the abundances of A, A, B, and B track the equilibrium distributions until the AA → φφ and BB → φφ annihilation processes freeze out at x 20 (= -9 -JCAP02(2016)052 x A ) and x 20 m A /m B (= x B ), respectively. 12 At that point, the A and B abundances rapidly decrease due to the asymmetry. Assuming m B m A as in eq. (1.3), then the effective symmetry imbalance temperature is set by x A , yielding T I m A /20. Switching to the y variable, the AA → Bφ semi-annihilation process freezes out around y F 20, and T I T F since x A > y F . Eventually, the B particles decay, transferring their asymmetry back to SM baryons. We are left with a relic DM abundance of A and a baryon asymmetry η b , which for this benchmark, match the observed values in our universe.
We now can understand why it is important to take T D > T I . If B were in chemical equilibrium with baryons all the way down to T I , then the asymmetry in A would be very suppressed. Taking eq. (3.2) with x A 20, we find η A 10 −18 , which is far too small to obtain a reasonable DM abundance. For this reason, we need to have B decoupling happen before a large A/A asymmetry develops. That said, even though B is chemically decoupled from SM baryons at T D , one can still think of the dark sector as being in quasi-equilibrium with baryons out to T F , since any B particles produced will eventually decay to SM baryons. In this way, the asymmetry in the dark sector effectively leaks into the baryon sector via the AA → Bφ process, and this leakage stops only when semi-annihilation freezes out.
Perhaps less obvious is why it is important to take T I T F , which requires m B m A since T I is set by the lower freezeout temperature between AA → φφ and BB → φφ. The reason is that if there were a large abundance of B at later times, then the AB → Aφ process would still be active during A freezeout, severely depleting the DM abundance. 13 We will show this effect numerically in figure 7(b) below.
In the extreme asymmetric limit and assuming T D T I ≥ T F , we can gain an analytic understanding of the DM and baryon abundances. In this limit, the final A abundance is negligible, so 1 2 Y A (∞) = η A (∞) at late times. After all remaining B particles have decayed, the observed baryon asymmetry today is Similarly, the present-day DM relic abundance is where s 0 = 2891 cm −3 is the entropy density today, and ρ C = 5.15 × 10 −6 GeV/cm 3 is the critical density. Thus, in this extreme asymmetric limit, finding the baryon asymmetry and the DM abundance reduces to finding Y A (∞). Moreover, in this limit, the value of Y A (∞) can be effectively determined by considering just the AA → Bφ process. As shown in appendix D, the relevant Boltzmann equation is: where λ = s/H, the Hubble scale H(T ) is evaluated at T = m A − 1 2 m B , and the equilibrium densities are defined in eq. (D.1). We used eq.  to the standard WIMP case [47], albeit with the replacement x → y, Y A stops following the equilibrium distribution after freeze-out, and the Boltzmann suppressed term approaches zero. Using the freeze-out approximation, the solution to eq. (3.12) is Note that the (Y eq A ) 2 /Y eq B term in eq. (3.12) scales like e −2y , compared to the WIMP case with (Y eq A ) 2 ∼ e −2x , which explains why freezeout happens at y F 20 instead of x F 20. We can gain further insights into Y A (∞) by considering limiting cases. In the limit of a large AA → Bφ cross section, the 1/Y A (y F ) factor drops out, and we recover the familiar WIMP approximation (3.14) Just as for standard WIMPs, y F has a logarithmic dependence on cross section [47]. In the limit of a small AA → Bφ cross section, the A and B particles are effectively decoupled after freeze-out, so the A abundance is set simply by Y A (y F ), or more accurately, the initial asymmetry at decoupling 1 2 Y A (∞) η A,D . For our benchmark scenario, we end up somewhere in between these two extremes, which is a hybrid of standard WIMP-like behavior (i.e. AA → Bφ freezeout) and asymmetric DM behavior (i.e. initial asymmetries).

Results for dark matter abundance
We now present numerical results for the DM abundance and compare to the analytic approximation in eq. 11.5 11.5 thermal freeze-out expectation (albeit with x → y compared to ordinary WIMPs). For very small values of σ AA→Bφ v , the A abundance saturates at the initial asymmetry.
In figure 7(b), we show how the DM abundance changes as a function of m B /m A . For m A m B < 2m A , Ω DM decreases slowly when decreasing the m B /m A ratio. Once m B < m A , there is a dramatic drop in Ω DM , which arises because the AB → Aφ process is still active and relevant for determining Y A (∞). It is important to note that when m B > m A , the equilibrium term (Y eq A ) 2 /Y eq B decreases sharply, faster than the other terms in the Boltzmann system, and therefore setting it to zero after y F is a valid approximation. This does not hold for m B < m A , and that is the reason for the difference in behavior between the two regimes. We are mainly interested in the extreme asymmetric limit where A and B decouple from the evolution of Y A (i.e. T I T F ), which is why we focus on the parameter space m A m B < 2m A . Again, for this mass range it is tempting to interpret B as a bound state of two As, though we will not restrict ourselves to that interpretation in this paper.

JCAP02(2016)052
In figure 8, we show slices through the {m A , m B /m A , σ AA→Bφ v , η tot } parameter space. The star marks the benchmark parameters in eq. (2.7). The red regions correspond to the observed DM abundance (Ω DM h 2 ∈ [0.115, 0.124]) and the blue regions correspond to the observed baryon asymmetry (η b ∈ [8.1, 9.5] × 10 −11 ). In figure 8(a), the desired DM abundance is achieved with a AA → Bφ cross section familiar from the standard WIMP case, with only a small variation with η tot . In figure 8(b), the required value of m A also has a weak dependence on η tot , due to the mild logarithmic dependence of y F on m A , just as in the WIMP case. Note that the resulting DM abundance is somewhat correlated with the total asymmetry, since the thermal cross section benchmark is at the transition between the asymmetric regime and the WIMP regime, as shown in figure 7(a).
Fixing η tot , we can highlight the mass and cross section dependence. In figure 8(c), we see the presence of two different regimes in agreement with eq. (3.13) and figure 7(a). In the low cross section limit, the DM abundance is dominated by the initial asymmetry, which is indicated by the vertical behavior at m A 15 GeV. At higher cross sections, m A has to rise linearly with the cross section in order to obtain the correct abundance using eq. (3.14). In figure 8(d), we also see two regimes which follow from figure 7(b). For m B /m A > 1, the DM abundance band has only a mild dependence on m A . As m B /m A → 1, the number density drops dramatically for fixed cross section, so the value of m A has to increase to compensate.

Fitting the Galactic Center excess
A number of studies have used the Fermi public data to identify a potential excess of gamma rays coming from the GC region [11][12][13][14][15][16][17][18][19][20][21][22]. The origin of this excess is as-of-yet unknown, but a tantalizing possibility is that this a signature of DM annihilation, though astrophysical explanations are also plausible [16,18,19,[65][66][67][68]. The DM interpretation is bolstered by the fact that the needed annihilation rate is consistent with that of a thermal WIMP, though there is recent evidence that the GC excess is better fit by a population of unresolved point sources [69,70]. In any case, typical asymmetric DM models do not predict this kind of indirect detection signal, though, unless there is a residual symmetric component (see e.g. [71,72]). Here, however, the semi-annihilation process AA → Bφ followed by B decaying to hadrons can give rise to an interesting gamma ray signal. In this way, the GC excess could be connected to the process of asymmetry sharing in the early universe.
The decays of B are prompt on cosmological time scales, so one can think of the AA → Bφ process as being a one-step cascade decay [51,52,55,73] (see also [74][75][76][77]) where the B subsequently decays to three quarks. Recall from eq. (1.4) that our benchmark decay modes are B → udd and B → cbb, though other flavor combinations are equally plausible. These quarks hadronize, and the resulting gamma ray spectrum comes primarily from neutral pions which decay via π 0 → γγ. The decays of φ are model dependent and may contribute to the gamma ray signal as well. For concreteness, we assume that φ dominantly decays to µ + µ − , 14 which only results in a small contribution to the gamma ray spectrum from final state radiation (FSR). For simplicity, our study ignores these possible FSR photons as well as photons from inverse Compton scattering or bremsstrahlung.

JCAP02(2016)052
To determine the gamma ray spectrum in the B rest frame, we use Pythia 8.185 [78] to construct a color singlet three quark final state, uniformly filling out the allowed three-body phase space. 15 We then use the default hadronization and decay model in Pythia to obtain the gamma ray spectrum dN/dE rest after letting all unstable hadrons decay. To go to the AA → Bφ rest frame, we boost the B particle by the gamma factor taking m φ = 0 for simplicity. The resulting gamma ray spectrum is given by (see, e.g., [43]) The produced flux of gamma rays as seen on earth is The J factor is the integral along the line of sight of the DM density. We adopt the normalization of ref. [22] where the region of interest (ROI) is |l| ≤ 20 • and 2 • ≤ |b| ≤ 20 • : The effective cross section accounts for the possibility that A comprises only a fraction of the total DM density: (4.5) Written this way, the resulting gamma ray spectrum depends on three parameters 6) and the choice of B decay channel.
To fit the GC excess, we use the procedure outlined in ref. [22]. Adopting the same notation as ref. [77], the chi-squared for a given parameter point is where C −1 ij is the inverse covariance matrix, obtained from ref. [22]. We show example fits (reasonably close to the best ones) of the photon spectrum in figure 9, for both the B → udd and B → cbb channels. Because the photons from heavy quark decays are softer than for light quarks, obtaining a good fit for the B → cbb channel requires a higher mass than for the B → udd channel. This is consistent with the observation for standard WIMP scenarios that b quark final states require higher DM masses than τ lepton final states [20]. Because we are considering parameter points for which m A m B , the effect of the boost is mild, pushing the peak of the (energy-squared-normalized) photon spectrum to slightly higher  values. Consistently with previous work, the best fitting (effective) cross section is close to the expected WIMP thermal cross section.
In figure 10, we show the parameter regions that give the best fit to the GC excess in the B → udd (purple) and B → cbb (orange) channels. In each plane, we first find the best fit χ 2 , and then show 1, 2, and 3 standard deviation contours for ∆χ 2 . In figure 10(a), we show the m A versus σ AA→Bφ v eff plane, leaving m B /m A fixed, showing that the best fit regions tend to have WIMP-like cross sections. In figure 10(b), we show the m A versus m B /m A plane, leaving σ AA→Bφ v eff fixed. The value of m B /m A determines the boost of B, and can be used to fine tune the gamma ray spectrum. Larger boosts (smaller mass ratios) are somewhat preferred by the fit.
In figure 10(a), we have superimposed the CMB limits discussed below in eq. (5.1). These limits are largely independent of m B and do not constrain the parameters in figure 10(b). As noted in refs. [3,79,80], models that tend to fit the GC excess are in some tension with the CMB limits, though there is still viable parameter space. The CMB limits are proportional to a parameter f eff that quantifies the fraction of energy that goes into electrons and photons. We compute this parameter below in eq. (5.3) and find f B eff = 0.17 if we only account for the B decay products, though we use the more conservative value f eff = 0.24 which assumes φ → µ + µ − .
In figure 11, we show the same parameter space as figure 10, but now letting the third parameter from eq. (4.6) float to give the best fit. In figure 11(a), we note that floating the ratio m B /m A extends the best fit to a wider range of values for both m A and σv eff but still within the vicinity of the benchmark in eq. (2.7). In figure 11(b), we see that the best fit for each operator is determined by a band in m A , and the change in the mass ratio and can be compensated somewhat by a change in the cross section.
Finally, returning to figure 2 from the introduction, we combine the analysis of the DM and baryon abundances in figure 8(c) and the GC best fit regions in figure 11(a). Note that JCAP02(2016)052 the abundances are given with respect to the actual cross section σ AA→Bφ v , while the GC fits and CMB bounds are with respect to the effective cross section σ AA→Bφ v eff . We made this hybrid choice to avoiding display a pathological region of phase space where one obtains a good fit to the GC excess with a small cross section but overabundant DM. As advertised, the benchmark parameters in eq. (2.7) yield a consistent cosmology and a plausible fit to the GC excess.

Additional constraints and signals
Having seen that we can achieve a viable asymmetric DM scenario with intriguing indirect detection signals, we briefly discuss possible additional constraints and signals.
• CMB heating bounds. The process AA → Bφ can occur in the early universe, even after thermal freeze-out. This residual semi-annihilation is constrained by limits on the power injected into the CMB through ionization [81]. The power in this case is parameterized by: where f eff is the efficiency factor and the effective cross section is defined in eq. (4.5). We can divide f eff into contributions from the B decay products and the (model-dependent) φ decay products.
Following the analysis of ref. [82], the efficiency factor from species X is: where we read off the values of f i eff , i ∈ {e ± , γ} from ref. [82] and estimate f p eff ≈ 0.2(f e + +e − eff + f γ eff ) following ref. [83]. For the B → cbb decay, we find f B eff = 0.17. For the φ decay, we assume that the dominant decay mode is µ + µ − , leading to f φ eff = 0.07. Current constraints from Planck [3] are shown in figures 10(a) and 11(a) above (for figures 10(b) and 11(b), the bounds are outside of the plotted region as they constrain m A 2 GeV). Note that the power injected depends directly on σ AA→Bφ v eff , and therefore directly impacts the GC excess best fit regions. For a fixed effective cross section, the CMB limits become less stringent as m A increases.
• Antiproton flux bounds. Though B has baryon number +1, we nevertheless expect to obtain antiprotons from the AA → Bφ process, since the B decay products will hadronize. 16 This additional antiproton flux can be tested in cosmic ray experiments like PAMELA [84] and AMS-02 [85]. The flux of antiprotons is given by [86] dφ p dK

JCAP02(2016)052
where K is the kinetic energy of the antiproton (a distinction important in the low energy limit), v p is the velocity of the antiproton, and R(K) is a best fit function that describes the propagation of the antiprotons throughout the galaxy (see [86]). As in the gamma ray case from section 4, we can extract the antiproton spectrum dN p /dK from B decays in Pythia and boost to the AA rest frame.
Various groups have derived antiproton bounds using AMS-02 data [87][88][89], typically showing results for WIMPs that annihilate to bottom quarks. To a reasonable approximation, the antiproton yield in B → cbb decays is comparable to that of a bottom quark, yielding around 0.3 antiprotons per decay. More accurately, a single B → cbb decay from AA → Bφ yields a factor of 2.5 fewer antiprotons than an equivalent energy bb pair. Thus, instead of evaluating eq. (5.4) directly, we can simply scale down the χχ → bb bounds by this factor. Taking the "Ein MED" bounds from [87], for m A = 60 GeV we estimate which is in some tension with the GC best fit region. That said, there is at least a factor of two or three astrophysical uncertainty in these bounds. In addition, while the CMB bounds are irreducible, in the sense that the same photons that explain the GC excess will inevitably correspond to power injected into the early universe, perturbing the CMB, the antiproton bounds are less directly tied to the gamma ray signal. For both of these reasons, we have opted not to show the antiproton flux bounds in our plots.
• Direct detection bounds. Because A has couplings to φ and φ couples to SM states, there will necessarily be a contribution to A-nucleus scattering from t-channel φ exchange. For specific models, such bounds might be relevant (and prospects for future direct detection experiments promising), though direct detection constraints can typically be avoided for two reasons. First, while φ needs to have large enough coupling to stay in thermal equilibrium with the SM, those required couplings are small from the perspective of A-nucleus scattering. Second, small increases in m φ can lead to large decreases in the A-nucleus cross section, since t-channel scattering typically scales like 1/m 4 φ at small recoil energies. For the specific model studied in appendix A, φ mixes with the SM Higgs after electroweak symmetry breaking with a mixing angle θ φh . The spin-independent scattering cross section of A through φ is (see, e.g. [57,90]) where λ A is the coupling of A to φ, v EW is the Higgs vacuum expectation value, f is a factor obtained from the different parton fractions which we take to be 0.35 [91], and m n is the nucleon mass. Taking m φ = 250 MeV as a benchmark, this cross section scales like where the baseline value saturates the LUX bound [92]. While this baseline mixing angle is small, φ still decays prior to BBN. Of course, larger mixing angles are allowed -18 -

JCAP02(2016)052
for larger m φ , though one should then account for φ decays in the GC excess analysis (see footnote 14).
• Flavor bounds. One aspect of this scenario that we have not delved into deeply is the generation of the B decay operators in eq. (2.4), and in particular their flavor structure. At the scale Λ (or below), some new heavy states are required to generate these operators, and those states could contribute to flavor-violating interactions. For any Bu c d c d clike operator and assuming that Λ-scale physics conserves CP, one expects that the strongest limits should come from meson mixing induced by those new heavy states [93][94][95]. Since all fields involved in the B decay are right-handed, stronger constraints from chirality-mixing operators can be avoided. Assuming a generic flavor structure for the new physics, the most constraining bound comes from the (c R γ µ u R ) 2 /Λ 2 operator, which is bounded by Λ ≥ 1200 TeV [94]. Taking m B = 60 GeV in eq. (2.5), this yield a B decay width of Γ B 1.2 × 10 −20 GeV. Interestingly, Γ B ∼ H occurs at a temperature of 121 MeV, above the beginning of BBN at T 10 MeV [46,61], so this flavor-safe scenario is indeed consistent cosmologically. Smaller B lifetimes (such as for the Λ = 300 TeV benchmark we use in eq. (2.5)) require some mild suppression of flavor violation among the heavy states, which is certainly plausible if Λ-scale physics is approximately flavor conserving.
• Collider searches for displaced jets. The operators that lead to B decay can also lead to B production at the colliders. For our B → cbb benchmark, this cross section is rather small at the LHC, but if the B → udd operator is present, then there is a more promising process: ud → Bd. 60 GeV and B decays to quarks, this yields a relatively low energy four jet final state, and it is questionable if such events could be seen over overwhelming QCD backgrounds. That said, plugging Λ = 30 TeV into eq. (2.5), the B has a lifetime of τ B 2.2 × 10 −11 sec, or a decay distance of cτ B 0.7 cm. Thus, the jets from B decays come from a (potentially very) displaced vertex, similar to the phenomenology of hidden valleys [96,97] and their variants (see e.g. [98][99][100]).
To get an estimate of the B production rate, the parton-level cross section in eq. (5.8) scales like To obtain the proton-proton cross section for the 14 TeV LHC, we integrate over the MSTW2008 LO parton distribution functions [101]: For the high luminosity LHC, with target luminosity of 3 ab −1 , we expect around 120 events. With a dedicated displaced jet trigger, one could hope to identify these events.

JCAP02(2016)052
This would be a distinctive signature for this scenario, especially if one could somehow verify that every B decay yields a baryon +1 final state. We leave a study of the LHC detection prospects to future work (see related studies in [102][103][104]).

Conclusions
In this paper, we presented a DM scenario where the interactions responsible for asymmetry sharing in the early universe are potentially visible today through indirect detection experiments.
The key novelty compared to other asymmetric DM scenarios is that the DM abundance is set by thermal freeze-out involving an unstable particle B rather than by the decoupling of high-scale interactions. Assuming that B chemically decouples above the DM mass, the parametrics of A freeze-out behaves much like a standard WIMP, albeit with a non-zero chemical potential and a rescaled freeze-out value x F → y F . Given this connection to WIMP physics, it is not surprising that we found a benchmark scenario with the right DM and baryon abundances. Intriguingly, this same benchmark yields a gamma spectrum from AA → Bφ semi-annihilation compatible with the GC gamma ray excess seen by Fermi. This work emphasizes that non-minimal asymmetric DM scenarios can produce interesting indirect detection signals without relying on a residual symmetric DM component. As shown already in refs. [41,43,44], asymmetric DM models with multiple stable states can yield indirect detection signals from semi-annihilation. Here, we emphasize the role that unstable dark sector states can play both in generating indirect detection signals as well as in sharing primordial asymmetries.
An interesting variant to our scenario is if B were stable, perhaps due to an additional Z 2 symmetry. In that case, both A and B would contribute to the DM abundance, with the relative ratio determined by the AA → Bφ process. As long as there is a sufficient abundance of A particles, then this semi-annihilation process could give rise to both a boosted DM signal from the final state B [105-108] as well as a standard indirect detection signal from the decays of φ. We leave a study of this scenario to future work.
Finally, the nature of DM cannot be determined through any single observation, and even if the GC excess is indeed due to DM (semi-)annihilation, one would want to determine the dark sector properties through a suite of other experiments. For this particular DM scenario, the hadronic B decays imply an irreducible cosmic ray signal from antiprotons and potentially antideuterons. More intriguingly, evidence for B particles might show up at the LHC or future colliders through displaced hadronic decays, yielding a visible portal to a rich dark sector.

A Concrete model
For simplicity in the text, we referred to A as being a single particle. Since the decay operators in eq. (2.4) imply that B must be fermion, though, A should be replaced by a fermion/boson system of ψ A /φ A . Motivated partly by supersymmetry but mainly by the resulting simplification of the Boltzmann equations, we take ψ A /ψ c A to be a Dirac fermion that is mass degenerate along with two complex scalars φ A /φ c A (i.e. two chiral multiplets with a holomorphic mass). We comment on the impact of splitting the ψ A /φ A masses below.
In the text, we also referred to φ as being a single particle and used the benchmark m φ = 250 MeV for studying CMB and direct detection bounds. Here we replace φ with two real fields, one scalar and one pseudoscalar. This can also be motivated by supersymmetry, but more relevant to our scenario, if φ were composed of a single real scalar field, then the annihilation ψ A ψ A → φφ would be p-wave suppressed. We consider these two fields to have somewhat heavier masses than considered in the text, in order to avoid additional bound state formation (see footnote 7). We will see that this has a negligible effect on the CMB bounds in eq. (5.3) and weakens slightly the direct detection bounds in eq. (5.6).
The field content for this model is shown in table 1, where B/B c is a Dirac fermion and φ is a complex scalar. We separate φ into its scalar component s and pseudoscalar component a, as This is similar to the axion portal [55,56], though we have suppressed a possible vacuum expectation value (vev) for φ in order to remain agnostic as to whether or not a is a pseudo-Goldstone boson from spontaneous symmetry breaking. The Lagrangian is then given by

JCAP02(2016)052
Note that we have introduced a mass splitting between the scalar and pseudoscalar modes, such that the decay s → aa is kinematically allowed. For simplicity, we have taken the Yukawa couplings to s and a to be the same, and assumed various relations among the scalar φ A and φ c A couplings. The potential term V (H, φ) includes terms that mix φ with the Higgs sector after electroweak symmetry breaking. Motivated by supersymmetry and the axion portal, we assume a two Higgs doublet model with H u,d , such that one possible mixing term is This interaction assures that both s and a can decay to SM states, though V (H, φ) generically contains cubic interactions such that the s → aa decay mode dominates. We ignore the induced vev of φ from this mixing, since it can be absorbed into a redefinition of the other couplings. Depending on the potential, the Higgs vevs can also contribute to the φ mass, though we expect this to be a small effect given the small mixing angle suggested by eq. (5.6). Without supersymmetry, some degree of fine tuning would be necessary to keep the s and a masses small given their large couplings to the A and B states. The baryon assignment of the different particles is set by Bu c d c d c /Λ 2 in eq. (A.4). So unlike in the SM, baryon number is not an accidental symmetry of this Lagrangian, though that conclusion might change depending on the precise dynamics present at Λ.
For the purposes of appendix D, the key feature of the coupling choices above is that the A scalars are treated symmetrically, such that any process involving φ A has a counterpart involving (φ c A ) † . This, along with the assumed mass degeneracy of the A fermions and A scalars, allows the full Boltzmann system to simplify to a two-particle system. We checked that all the relevant DM interactions from this Lagrangian include an s-wave term (see e.g. [109,110]), which is necessary to achieve the desired cosmology. For example, an s-wave annihilation channel for the A fermions is possible via ψ A ψ A → sa.
Introducing mass (or coupling) splittings between ψ A and φ A would add new parameters to the model studied and complicate the Boltzmann equations. The dominant effect of such a splitting is to change the resulting populations of the ψ A and φ A components. To maintain roughly the same phenomenology as presented in the text, we need m ψ A /m φ A to be O(1), such that the lighter DM component is still Boltzmann suppressed at freezeout. We also need to satisfy the mass ordering This first inequality comes from the requirement of not depleting A particles by scattering with B particles prior to freezeout (see footnote 13), and the second inequality ensures that the semi-annihilation process ψ A φ A → Bφ is kinematically allowed even at threshold. In addition to the possibility that B might correspond to a bound state of A particles, the φ mediator might give rise to additional bound states of A and/or B. For simplicity of our analysis, we want to avoid the possibility of φ-mediated bound states, which would further complicate the Boltzmann analysis. Note that only the scalar s mediates a 1/r Yukawa potential between the A particles, whereas the pseudoscalar a can only mediate a 1/r 3 spindependent potential among the A fermions. For our benchmark, with O(1) φ couplings and DM masses around 50 GeV, we can avoid bound state formation if m s O(1 GeV).
The presence of both s and a turns out to have relatively little effect on the phenomenology presented in the body of the paper. There are now two semi-annihilation processes relevant for indirect detection, ψ A φ A → Ba and ψ A φ A → Bs. As long as 2m µ < m a < 3m π -22 -

JCAP02(2016)052
(such that a dominantly decays as a → µ + µ − ) and m s > 2m a (such that s dominantly decays as s → aa → µ + µ − µ + µ − ), then the resulting photons from muon FSR give subdominant contributions to the photon flux as needed for section 4. The CMB bounds in eq. (5.3) only have a mild dependence on the injected muon spectrum so are largely unchanged. Raising m s weakens the spin-independent direct detection bound in eq. (5.6) and spin-dependent bounds from a exchange are typically subdominant.
Finally, we note that similar phenomenology could be achieved by replacing φ with a U(1) gauge field that kinetically mixes with SM hypercharge [48][49][50]. The main challenge for using this hypercharge portal is that the B decay operator in eq. (A.4) would then require an insertion of the U(1) -breaking Higgs field, making it an even higher dimension operator.

B Chemical potential analysis
To determine the abundance of DM at the decoupling temperature T D , we have to consider the chemical potentials of all relevant SM and dark species. For simplicity, we work in a regime where the sphaleron process becomes inactive prior to T D , such that we can treat the baryon and lepton asymmetries as being independently conserved during B decoupling. The states in the dark sector then carry effective baryon number consistent with the operators in eq. (2.4).
In the early universe, SM interactions guarantee chemical equilibrium among SM particles. After electroweak symmetry breaking, the chemical potential relations can be written as [62,111] where µ u , µ d , µ ν , and µ e refer to the chemical potentials of each flavor of up-type quark (u L and u R ), down-type quark (d L and d R ), left-handed neutrino (ν L ), and charged lepton (e L and e R ). Note that eq. (B.1) is enforced by W ± exchanges (with chemical potential µ W ) and eq. (B.2) is imposed by the sphaleron process. At temperatures below the top quark mass, the neutrality condition imposes The above relations imply which allows us to write the chemical potential for SM baryons as For a species i with g i degrees of freedom and mass m i , the relation between the charge density and the chemical potential µ i at temperature T is [112] (see also [40]) where we are assuming µ i T with (B.10) The function f (x) takes into account the Boltzmann suppression of particle i. In the limit Using eqs. (B.5)-(B.10), we can write the asymmetries as η i = X i (n i − n i )/s, where X i denotes the baryon charge of species i. This leads to where g q is the number of degrees of freedom for each quark flavor, i.e. g q = 3 from color (since the left-and right-helicities are already included in eq. (B.5)). For the system described in appendix A, The relations above are valid for temperatures T ≥ T D , where T D is the B decoupling temperature. Using the total asymmetry η tot in eq. (3.1) can be written as which is the basis for eqs. (3.2)-(3.4). From eq. (B.9) and the relation above, it is also possible to write the DM chemical potentials as functions of η tot as where we used s = 2π 2 g * s T 3 /45. The process Bq → qq, where q is a quark consistent with the operators in eq. (2.4), is efficient at high temperatures, and must be included when determining the B decoupling temperature T D . From the amplitude expressed in terms of the Weyl wavefunctions x i and y i , we find where s is the square of the center-of-mass energy, and we assume m B m q . The differential cross section with respect to the Mandelstam t variable is where to simplify the phase space factor we use the function Integrating the cross section over t, we find σ = s 64πΛ 4 . (C.5) At sufficiently high temperatures, we can ignore the chemical potentials of B and the quarks. For a generic scattering process 12 → 34, the Møller thermal cross section is [113] σv Mol 12→34 = 1 n eq 1 n eq 2 g 1 d 3 p 1 (2π) 3 where the Møller velocity is defined such that the product v Mol n 1 n 2 is Lorentz invariant, and In our case, n eq 1 = n eq B0 is the thermal Boltzmann distribution for B, where K i is the modified Bessel function of order i, and n 2 = n eq q0 is the relativistic distribution for quarks, where here we have corrected for the Fermi statistics of the relativistic quarks. The relative velocity v Mol is obtained from p 1 · p 2 = v Mol E 1 E 2 . The integral in eq. (C.6) can be reparameterized following the annihilation example of ref. [113] as

D Boltzmann equation details
In going from the full model in appendix A to the simplified discussion in section 2, we asserted that the four-particle system X ∈ {ψ A , φ A , φ c A , B} could be reduced to an effective two-particle system. Indeed, this is possible if all states type A are mass degenerate and have equal chemical potentials such that µ ψ A = µ φ A = µ φ c † A ≡ µ A . We now derive the full Boltzmann system and show how this reduction occurs.
Unlike the discussion in section 3.3 where we introduced the y variable, here we stick to the notation x = m A /T . The equilibrium distributions for a state X are given by Y eq X (x) = Y eq X0 (x) exp(µ X /T ), (D.1) where µ X is the chemical potential of X and we assume µ X T which is a valid assumption until freezeout. In the non-relativistic limit (m X T ), bosons and fermions follow the same distribution [47], and the equilibrium functions denoted with a subscript "0" are Y eq X0 = g X g * s 45 4π 4 x 2 K 2 (x) . (D.2) We use the notation λ = s/H(m A ), where s is the entropy of the universe, H the Hubble constant at the temperature T = m A , and g * s is the effective number of relativistic degrees of freedom at the same temperature [47]. The chemical potentials satisfy 2µ A = µ B when the semi-annihilation process AA → Bφ is in equilibrium. We also have µ A = −µ A and µ B = −µ B when the annihilation processes AA → φφ and BB → φφ are in equilibrium.
Here we assume µ φ = 0, since the φ is in thermal equilibrium with the SM and there is no symmetry forbidding interactions like φφ → φφφ.
By assumption in appendix A, the couplings involving φ A and φ c † A are identical, so we know that at all temperatures T . Therefore, processes that involve φ A , such as σ ψ A B→φ A φ v , will simply get a factor of 2 from including the equivalent process with φ c † A . With this simplification, the -26 -

JCAP02(2016)052
Boltzmann equations for the {ψ A , B} system are: The Boltzmann equations for φ A can be easily extrapolated from those of ψ A . Here, we have neglected the conversion process AA → BB since annihilation dominates over conversions for determining the freezeout dynamics of near-mass states (recall that our mass range of interest is m A m B < 2m A ). At high enough temperatures (though not so high that boson/fermion statistics matter), the equilibrium distributions satisfy where the factor of 2 comes from the number of degrees of freedom. In general, though, the evolution of Y ψ A and 2 Y φ A will not remain identical, except for special choices of the couplings in eq. (A.3). To see which combination is needed, consider the evolution of -27 -

JCAP02(2016)052
In addition to the factors of 2 coming from including φ A ↔ φ c † A , there is an additional factor of 2 in the second line coming from averaging over polarizations in the cross section, since There is a factor of 1/4 in the second line from the ratios of the equilibrium distributions since Y eq ψ A = 2Y eq φ A . In the third line, the factor of 1/2 comes from averaging over spins in the cross sections, and the factor of 4 comes from the ratio of equilibrium densities. From the semi-annihilation terms alone (i.e. second and third lines), it is now clear that if Y eq ψ A = 2Y eq φ A at early times, then this relation will continue to hold at later times. For the annihilation terms (i.e. first line) to maintain the same relation, we simply need to choose couplings such that This in turn sets the needed couplings in eq. (A.3). Of course, for more general couplings, one can simply solve the full set of Boltzmann equations. Assuming eq. (D.9) holds and plugging the relation Y ψ A = 2Y φ A into the Boltzmann equations, we arrive at an effective two-particle system with {A, B}: with conjugate equations for A and B. Here we refer to Y ψ A as simply Y A for simplicity, but we maintain old notation for the cross section in order to maintain the correct factors of 2 in averaging over polarizations. As one can check, the B decay term and the Bq → qq scattering term depend explicitly on the chemical potentials, as expected from the discussion in section 3.2. Below T D , though, when the scattering term can be neglected, the chemical potential factors cancel in the Boltzmann system, since we take ratios of the equilibrium densities except for the decay -28 -

JCAP02(2016)052
term. In other words, after T D , we could make the replacement Y eq X → Y eq X0 , as defined in eq. (D.2). Numerically this is equivalent to setting the chemical potentials solely as initial conditions, as given by, say, eq. (3.8).