Indirect detection of dark matter with (pseudo)-scalar interactions

Indirect detection is one of the most powerful methods to search for annihilating dark matter. In this work, we investigate the impact of non-perturbative effects in the indirect detection of dark matter. For this purpose we utilize a minimal model consisting of a fermionic dark matter candidate in the TeV mass range that interacts via scalar- and pseudo-scalar interactions with a massive scalar mediator mixing with the Higgs. The scalar interaction induces an attractive Yukawa potential between dark matter particles, such that annihilations are Sommerfeld enhanced, and bound states can form. These non-perturbative effects are systematically dealt with (potential) non-relativistic effective field theories and we derive the relevant cross sections for dark matter. We discuss their impact on the relic density and indirect detection. Annihilations in dwarf galaxies and the Galactic Center require special care and we derive generalized $J$-factors for these objects that account for the non-trivial velocity dependence of the cross sections in our model. We use limits on the gamma-ray flux based on Fermi-LAT observations and limits on the rate of exotic energy injection from Planck to derive bounds on the parameter space of the model. Finally, we estimate the impact that future limits from the Cherenkov Telescope Array are expected to have on the model.


Introduction
The advent of precision cosmology, mainly driven by increasingly accurate observations of the Cosmic Microwave Background (CMB), has made stunning insights into the structure and the dynamics of our universe at the largest scales possible.Among the most startling results of this program is the realization that the energy content of the universe is heavily dominated by two components that are only poorly understood: dark matter (DM) and dark energy.The concordance cosmological model, known as ΛCDM, requires that ∼27% of the energy budget of the universe is in the form of cold DM to explain the CMB anisotropies measured by the Planck satellite [1].This constitutes a serious challenge for particle physics since none of the known particles can account for it.
Both on the side of ultraviolet complete theories and simplified models, a plethora of particle physics models have been proposed that contain viable dark matter candidates (see e. g. [2,3]): stable, or very long-lived, particles that are electromagnetically neutral and interact only (very) weakly with the visible sector.The possibility to observe DM particles with particle-physics-inspired strategies requires additional interactions between the visible and the dark sector beyond gravity.One organizing principle that allows to classify different DM candidates, is the production mechanism in the early universe.This is particularly helpful if one is interested in the phenomenology of the DM particles since different production mechanisms require vastly different masses and couplings.One of the classic solutions is the thermal production of DM via freeze-out from the SM bath [4].This has inspired a wide range of experiments, such as direct detection, indirect detection, and collider searches, see e. g. [5][6][7][8][9] for some recent reviews.Out of these approaches, indirect detection stands out as being the most immediately connected to the production process.It aims to observe the continued annihilation of DM particles in the late universe after freeze-out and picks out this signal from backgrounds of astrophysical origin.There is a vast range of indirect searches that target very different messengers such as gamma rays [10][11][12][13][14][15], antiprotons [16][17][18], positrons [19][20][21] and neutrinos [22][23][24].All of these have their individual strength and limitations.The comparative ease of detection, and the fact that photons travel unaffected from the source to Earth, make gamma rays particularly interesting.With the Cherenkov Telescope Array (CTA) a next-generation gamma-ray telescope is on the way [25].Therefore, it is timely to reassess the impact of existing indirect detection limits on dark matter produced by freeze-out and compare this with the potential of CTA observations.A particularly interesting scenario for indirect detection is so-called secluded DM [26,27], i. e. DM annihilating into states beyond the SM.This naturally allows to decouple the relic density from other observables and thus avoids stringent bounds from direct detection experiments [28][29][30] and the LHC [31][32][33][34].When these new states are light compared to the DM mass, special care is needed for the computation of the annihilation rates since they can lead to long-range forces between the annihilating DM particles.The presence of such a long-range interaction leads to a modification of the annihilation rate by the Sommerfeld effect [35][36][37] and, if the interaction is attractive and strong enough, it allows for the existence of bound states [38,39].
On the one hand, the inclusion of Sommerfeld factors has become over the last years a somewhat standard ingredient in both the extraction of the DM relic density and the indirect detection limits [36,[40][41][42][43][44].On the other hand, accounting for bound-state formation is less established and the methods are still being developed.In this case, a large body of results focuses on one aspect at a time (DM energy density [39,[45][46][47][48][49][50][51][52][53][54][55][56][57][58][59] or experimental signatures [38,[60][61][62][63][64] respectively), whereas there are a few studies in the literature that consistently combine the impact of such non-perturbative effects on the cross sections for the thermal freeze-out and the indirect detection.Relevant examples are found in refs.[50,[65][66][67], where the impact of a vector mediator, i. e. a dark photon which undergoes kinetic mixing with the U(1) Y of the SM, is scrutinized.In this work, we aim to fill the gap for a class of DM models that feature self-interaction as mediated by a scalar particle.More specifically, we consider both scalar and pseudo-scalar interactions with a Dirac fermion dark matter.Using non-relativistic effective field theory, we derive the cross sections that fix the DM energy density and play a role in DM annihilations in astrophysical environments, such as the Galactic Center (GC) and dwarf spheroidal galaxies (dSphs).These are then employed to map out the cosmologically viable parameter space for the model and combined with present and prospective limits from the gamma-ray telescopes Fermi-LAT and CTA, respectively.
The paper is structured as follows.In section 2 we describe the DM model and its relevant energy scales.Then, the cross sections and decay widths are presented in section 3 within the framework of non-relativistic effective field theories.The Boltzmann equation that we use to extract the DM energy density, together with a few highlights about the impact of Sommerfeld and bound-state formation, are given in section 4. The scrutiny of various astrophysical and cosmological environments that can constrain present-day DM annihilation is addressed in section 5. Here, the current bounds and future prospects for the model, which include the next-generation instrument CTA, are presented.Conclusions and outlook are offered in section 8.

DM model and energy scales
In this section, we describe the field content of the DM model and discuss the relevant energy scales.We assume DM to be a Dirac fermion that carries no charge under the SM gauge group, namely, it is an SM singlet.However, we introduce an interaction between DM fermions as mediated by a scalar particle via Yukawa-type interactions.The scalar mediator of the dark sector couples to the SM via a Higgs portal interaction.The Lagrangian density of the model reads where X is the DM Dirac field and ϕ is a real scalar.The scalar self-coupling is denoted with λ ϕ , whereas the scalar and pseudo-scalar couplings with the fermion are g and g 5 respectively.For simplicity, we assume the scalar self-coupling to be negligible such that plays no role in our analyses. 1The mass of the scalar mediator m ϕ is assumed to be smaller than the DM mass M , in order to allow long-range effects.In this work, we restrict to the case where the scalar coupling is larger than the pseudo-scalar coupling, namely α ≡ g 2 /(4π) ≫ α 5 ≡ g 2 5 /(4π).In doing so, we ensure that the dominant nonperturbative effects originate from a scalar-type interaction, that induces an attractive potential, and we can largely neglect the mixed scalar-pseudo-scalar and pure-pseudo-scalar induced contributions [56,68].
For the main scope of this work, namely exploring non-perturbative effects in indirect detection of dark matter, it is sufficient to work within a simplified model realization with a scalar particle mediating interactions between DM particles [69,70].The mass parameters and couplings of the Lagrangian given in eq.(2.1) are taken as free and we do not address here the fermion and scalar mass generation of the dark sector (see e. g. [71,72] for a simplified model with two mediators, a vector and a scalar, where the gauge invariance and spontaneous symmetry breaking in the dark sector is fully accounted for).It is worth mentioning that the mediator sector can be extended to have a richer set of interactions [68,73,74].For example, an interaction of the form ρ ϕ ϕ 3 with a dimensional coupling could exist, which enables additional bound state formation processes [54,74].
The dark and visible sectors are connected via Higgs-portal interactions, which read as follows where H is the SM Higgs doublet.We have shifted the scalar mediator field such that it has no vacuum expectation value.In the following, we shall take the quartic coupling λ ϕh to be negligible and only retain the dimensional coupling to induce the mixing between ϕ and the Higgs boson. 2 We decided to induce the mixing between the scalar and the Higgs through a dimensional coupling µ ϕh (and with v ϕ = 0) in order to keep the number of relevant free parameters at a minimum.Note that we could have also set µ ϕh = 0 and kept the quartic mixing λ ϕh alongside the quartic dark-scalar coupling λ ϕ .This choice would introduce a vev for the scalar v ϕ ̸ = 0 in analogy to the Standard Model Higgs sector.The mixing term in the Lagrangian after EWSB as well as the definition of the mixing angle δ remain the same if one substitutes µ ϕh → v ϕ λ ϕh .We have checked that the phenomenology of both scenarios, together with the resulting dark scalar and Higgs boson masses, is equivalent up to O sin 2 δ of the scalar mixing angle, which we will neglect in the following due to its smallness.The scalar mixing and the implications for the indirect detection are addressed in section 5.3.1.Fermion-antifermion pair annihilation into dark scalars is the key process both for fixing the DM energy density via the freeze-out mechanism and for providing the astrophysical signatures of indirect detection.X X annihilations in the early universe, as well as laterstage annihilations, either during the formation of the CMB or today, occur in the nonrelativistic regime.Here fermions and antifermions move with relative velocities v rel ≪ 1.In this kinematical region, also referred to as near-threshold, repeated scalar exchange between DM particles with v rel ≲ α induce ladder diagrams that are not suppressed and need to be resummed.This resummation manifests itself in the dynamical generation of bound-state poles of order M α 2 at negative energies (below threshold) and a continuum spectrum of scattering states at positive energies (above threshold).The typical momentum exchanged between the fermion-antifermion pairs for v rel ∼ α is M α, which is of the order of the inverse Bohr radius of the bound states.The dark matter mass and the dynamically generated scales are the more separated the smaller α is, namely M ≫ M α ≫ M α 2 .In the following, we refer to these energy scales as hard, soft, and ultrasoft respectively.Another relevant energy scale is the mediator mass m ϕ .We will take it to be smaller than the DM mass M and at most as large as the soft scale M α.Thus, long-range interactions typically still occur, however, there are sharp deviations from the Coulombic regime that manifest in various ways.From a phenomenological point of view, we are mainly interested in m ϕ ≥ 2m π 0 .Lighter mediators do not produce an appreciable amount of γ-rays, which is the main observation signature considered in this work.However, lighter mediators can have implications for other observables such as the positron flux or the DM mass distribution in dense halos, see e. g. [68,75].The main difference to massless mediators is that DM self-interactions are screened when the soft scale M α is of the order of the mediator mass m ϕ .For scattering states subject to an attractive potential, the Sommerfeld enhancement is typically reduced and does not monotonically grow at small velocities.Moreover, a rich resonance structure is found that depends on the mediator and dark matter mass ratio m ϕ /M [35][36][37].As far as bound states are concerned, there is even a condition that has to be satisfied for their existence of at least one bound state, namely 1.6 m ϕ < αM [76].In order to make contact with previous work on the subject, we define the two dimensionless variables for later use The latter variable involves the Bohr momentum M α/2 or, equivalently, the inverse Coulomb Bohr radius a 0 ≡ 2/M α.The Coulombic limit is recovered for large velocities with respect to α, namely ζ ≪ 1, and for mediator masses much smaller than the soft scale M α, i. e. ξ ≫ 1.
Finally, thermal scales are relevant for the freeze-out of DM particles in the early universe.Most notably, there is the temperature of the thermal plasma T , which at the thermal freeze-out of massive relics is smaller than the DM mass, T ≪ M .However, the temperature can still be of the order of the dynamically generated soft and ultrasoft scales M α and M α 2 .In this work, we include thermal effects due to the medium assuming that the temperature is about the ultrasoft scale M α 2 or smaller.This implies that thermal effects do not enter the potential, which may be taken as an in-vacuum Yukawa potential.For thermalized dark matter fermions, one finds v rel ∼ T /M , so that their typical momentum is M v rel ∼ √ M T ≲ M α, which implies v rel ≲ α.For the temperature range considered in this work, it also holds that √ M T ≳ M α 2 .These conditions qualify M v rel as a soft scale and M v 2 rel ∼ T as an ultrasoft scale.Therefore, the hierarchy of energy scales, that we show in figure 1, is The mediator mass is taken to be m ϕ ≲ M α.Because of the assumption of the scalar sector interactions, namely no trilinear interactions and λ ϕ ≪ 1, thermal masses for the mediator are not considered in the hierarchy of scales (the condition T ≲ M α 2 would make it small in any case).In the following section 3, we discuss the non-relativistic EFTs that we use to compute near-threshold observables.Indeed, the use of the Lagrangian eq.(2.1) is in general unpractical since all energy scales are entangled in the amplitudes.

Non-perturbative effects in NREFTs
In the following, we find it convenient to express pairs in a scattering state with (X X) p , where p = M v rel /2 denotes the momentum of the relative motion, whereas a fermionantifermion pair in a bound state is indicated with (X X) n .In order to shorten the notation, n ≡ |nlm⟩ stands for the whole set of quantum numbers of a given bound state.In this section, we summarise the non-relativistic effective field theories that are useful to compute the main observables of interest.These include DM pair annihilations into light scalar mediators, which can occur both for above-threshold scattering states (X X) p and bound states (X X) n , and bound-state formation.We shall compute the relevant cross sections, as well as decays widths, in the framework of a potential non-relativistic EFT, dubbed as pNRYγ 5 , which is obtained from the model Lagrangian (2.1) via a two-step matching procedure [56,77,78].First, the modes with energy/momentum of the order of the hard scale M are integrated out, which results in NRY γ 5 .In this theory, the soft and the ultrasoft scales are still intertwined.In the second stage, the typical scale of soft momentum exchanges between non-relativistic fermions and antifermions is integrated out, together with the mediator mass m ϕ , that results in the lowest-lying EFT, namely pNRY γ 5 .
Here, the degrees of freedom are fermion-antifermion pairs that interact only with ultrasoft scalar mediators and obey a Schrödinger equation with a Yukawa potential.The towers of EFTs are shown in figure 1.
Here we discuss the non-relativistic EFT framework and its main ingredients for a self-contained discussion.The details on the derivation of the EFTs have been addressed in refs.[56,[77][78][79] and more recently in ref. [80] for the inclusion of thermal scales.
3.1 NRY γ 5 and pNRY γ 5 NRY γ 5 describes non-relativistic fermions and antifermions that interact with a scalar.The set of interaction vertices is comprised of the bilinear sector of the corresponding effective Lagrangian.Moreover, fermion-antifermion annihilations into light mediators are accounted for by four-fermion operators [81].The NRY γ 5 Lagrangian reads schematically where ψ is the two-component Pauli spinor that annihilates a dark matter fermion, whereas χ † is the Pauli spinor that annihilates an antifermion.Here, scalar particles have energy and momenta smaller than M .The fermion bilinear at order 1/M with the matching coefficients set to their tree-level values reads [56,79] for the non-relativistic fermion field, whereas for the antifermion.We notice that the leading fermion-pseudo-scalar interaction is suppressed with respect to the leading fermion-scalar interaction by k/M , where k is the soft or ultra-soft momentum carried by the field ϕ.The operators proportional to one power of g 5 are parity-violating, and the notation [∇ϕ] stands for the derivative acting on the scalar field only.The four-fermion Lagrangian L 4-fermions in eq.(3.1) is especially relevant because it accounts for DM fermion-antifermion annihilations into light scalars.Through the optical theorem, pair annihilations correspond to the imaginary parts of four-point Green's functions and, at leading order in the couplings, describe the 2 → 2 process X X → ϕϕ.We adopt the spectroscopy notation from NRQED/NRQCD [81,82] so that one can classify the annihilations in terms of the total spin S of the pair, the relative angular momentum L and the total angular momentum J.The operators and matching coefficients of the fourfermion Lagrangian are then labeled with 2S+1 L J (cf. eqs.(3.4) and (3.6)).This will turn out useful for assigning the corresponding Sommerfeld factors.In the following, we list the matching coefficients of the four-fermion operators by retaining the full dependence on the mass ratio m ϕ /M , and hence we generalize the results presented in ref. [56].In our setting of the hierarchy of scales, such corrections are modest (the largest value for the mediator mass is taken to be the soft scale M α).The general result with the m ϕ dependence could turn out to be useful if one aims to consider larger mediator masses for this model. 3n the following, we explicitly list the full set of four-fermion operators that can be written up to order 1/M 4 .The corresponding matching coefficients have to be determined and depend on the DM model.There are two independent dimension-6 operators that can account for velocity-independent contributions to DM annihilation, and they read [81] (L 4-fermions where f ( 1 S 0 ) and f ( 3 S 1 ) are the matching coefficients of the spin singlet and spin triplet operators Dimension-8 operators comprise derivatives acting on the DM fermion and antifermion fields and hence have a dependence on the DM velocity.Schematically one may write the corresponding Lagrangian in terms of operators, which are indicated with O( 2S+1 L J ) and P( 2S+1 L J ), and the corresponding matching coefficients f ( 2S+1 L J ) and g( 2S+1 L J ), respectively, as follows [81] (L 4-fermions where the operators are ) Here ↔ ∇ is the difference between the derivative acting on the spinor to the right and on the spinor to the left, namely χ † ↔ ∇ψ ≡ χ † (∇ψ) − (∇χ) † ψ.The notation T (ij) for a rank 2 tensor stands for its traceless symmetric components, The non-vanishing matching coefficients are found to be ) ) where the auxiliary functions F(r) and G(r) with r ≡ m ϕ /M are The vanishing imaginary parts for f ( 3 S 1 ) and f ( 1 P 1 ) can be understood in terms of symmetry arguments [56].The pseudo-scalar interaction in eq.(2.1) violates parity, however it still preserves the charge conjugation symmetry.Accordingly, only some combinations of the spin and angular momentum of the annihilating fermion-antifermion pairs are allowed. 4o obtain NRY γ 5 in eq.(3.1), we have integrated out the hard scale M , which is reflected by 1/M powers appearing in the effective Lagrangian.The next scale one has to integrate out according to the scale hierarchy in eq.(2.4) is the soft scale M α, or the typical inverse distance between X and X in a pair.In this work, we allow the mass of the scalar mediator to be as large as the soft scale M α, which also corresponds to the momentum transfer of potential exchanges among fermion-antifermion pairs (this is, in general, true for bound states, whereas it depends on the value of v rel for the scattering states).Hence, m ϕ is integrated out at the same time with the soft scale M α which results in a Yukawa potential experienced by X X pairs.In the situation where m ϕ ≪ M α one naturally recovers a Coulomb-like regime and the corresponding potential.As mentioned earlier, we can set the thermal scales to zero in the matching between NRY γ 5 and pNRY γ 5 because of the assumption T ≲ M α 2 .
To obtain pNRY γ 5 one implements a projection of the NRY γ 5 Hamiltonian onto the particle-antiparticle sector via where i, j are Pauli spinor indices, while the state |ϕ US ⟩ contains an arbitrary number of scalars with energies much smaller than M α, including thermal scales, and no heavy fermions/antifermions.Here φ ij (t, x 1 , x 2 ) is a bilocal wave function field representing the X X system.The Lagrangian is then written in terms of the relative distance of the pair r = x 1 − x 2 and its center-of-mass coordinate R = (x 1 + x 2 )/2 [56,78] where the square brackets in the second line of eq.(3.20) indicate that the spatial derivatives act on the scalar field only, which is multipole expanded.To avoid cluttering the notation, we suppress the spin indices of the bilocal fields that are contracted with each other.Each term in the Lagrangian (3.20) has a well-defined scaling.The time derivative scales as ∂ 0 ∼ M α 2 ∼ T , the inverse relative distance and the corresponding derivative obey r −1 , ∇ r ∼ M α (the mediator mass can as well define the soft scale whenever m ϕ ∼ M α), whereas the scalar field and the center-of-mass derivative satisfy gϕ, ∇ R ∼ M α 2 ∼ T .Indeed, the active dynamical scales in the so-obtained EFT are the ultrasoft scale and the temperature.
The potential is understood as a matching coefficient of pNRY γ 5 and it emerges as such because the typical distance between fermion-antifermion pairs is removed from NRY γ 5 .In general, the potential comprises both a real and an imaginary part, where the latter is related to physical processes X X → ϕϕ (see section 3.2).The potential is organized as an expansion in the coupling α and 1/M with the static Yukawa potential The contributions to the potential from the mixed scalar-pseudo-scalar and pure pseudoscalar diagrams have been scrutinised in ref. [56], and they belong to the 1/M -suppressed terms V (1) and V (2) , respectively.Using power counting of the low-energy theory, one can show that these terms scale as M α 3 (g 5 /g) and M α 4 (g 5 /g) 2 compared to the leading Yukawa potential in eq.(3.22), which scales as M α 2 .In this study, we work at (moderately) weak coupling, α < 0.25, and consider g 5 < g.Therefore, we can neglect pseudo-scalar effects in the potential to a good approximation.The same assumptions also guarantee that the corresponding pseudo-scalar-induced ultrasoft transitions among pairs are suppressed with respect to those displayed in the second line of eq.(3.20).Pure-scalar coupling diagrams contribute to V (2) and are α 2 -suppressed with respect to the leading Yukawa potential.Finally, because of T ≲ M α 2 , the potential does not get thermal contributions at any order.In what follows, we retain only the leading Yukawa potential (3.22) when solving the Schrödinger equation to find the corresponding wave functions of scattering and bound states.
The ultrasoft interactions in the second line of eq.(3.20) differ qualitatively from the potential-NREFT of vector mediators, namely pNRQED and pNRQCD [83].One finds a non-vanishing monopole term at order r 0 in the multiple expansion, whereas the electricdipole vertex, which is linear in the relative coordinate r, is absent.To capture the leading non-trivial dynamics among scattering and bound-states, one has to consider terms at order r 2 and find the quadrupole contributions that scale as M α 4 .The derivative interaction with ∇ r in the second line of eq.(3.20) arises from the 1/M 2 spin-independent operator of NRY γ 5 , and it is of the same order as the quadrupole term. 5

Heavy pair annihilations
Let us first discuss the annihilation of heavy pairs.This process is responsible for the depletion of the DM particles into lighter scalar mediators.Annihilations can equally occur for a particle-antiparticle pair in a scattering state or a bound state.At leading order in the coupling the processes are (X X) p → ϕϕ and (X X) n → ϕϕ, and are described by the same set of operators in V (2) , the only difference being on the two-particle state one projects them onto.
In NRY γ 5 , heavy pair annihilations are accounted for by local four-fermion operators, cf.eqs.(3.4) and (3.6).Accordingly, in pNRY γ 5 , the four-fermion operators generate local terms, which are proportional to a delta function δ 3 (r) in the imaginary part of the potential V (p, r, σ 1 , σ 2 ).We write the annihilation term from the Lagrangian density of pNRY γ 5 as follows [56,84,85] where S is the total spin of the pair, while T ij SJ and Ω ij SJ are spin projector operators (cf.e. g. [84,85]).The matching coefficients are given in eqs.(3.14)-(3.17).
To compute the annihilation cross section we use the optical theorem and consider the scattering amplitude from an initial to a final scattering state |p, 0⟩, where we have set the center of mass momentum to zero.By inserting the vertices that appear in eq.(3.23) in the self-energy of a scattering state, the spin-averaged annihilation cross section reads The dimension-8 derivative operators have produced two powers of momentum, that together with two inverse powers of the DM mass, gives p 2 /M 2 = v 2 rel /4 in eq.(3.24).The Sommerfeld factors S l ann (ζ, ξ) correspond to the wave function of the fermion-antifermion pairs in a scattering state evaluated at the origin, i. e. S 0 ann , where R l is the radial wave function with angular momentum l and p = M v rel /2.We stress here that, by computing the annihilation cross section for scattering states in pNRY γ 5 , the resummation of multiple soft-scalar exchanges (ladder diagrams) is already taken care of, since pNRY γ 5 is a quantum field theory for interacting pairs.The factorization of the different energy scales thus manifests in eq.(3.24): the hard dynamics is contained in the matching coefficients of NRY γ 5 , whereas the soft dynamics is contained in the wave-function squared, namely the Sommerfeld factors.With the exception of the Coulomb limit, m ϕ → 0, the Sommerfeld factors have to be numerically evaluated.We follow the procedure outlined in [87] and [45].
The complementary manifestation of multiple soft scalar exchanges between nonrelativistic fermions and antifermions is the appearance of bound states below threshold.Differently from the Coulomb case, there is a condition on the existence of bound states that depends on the relative sizes of the would-be Coulombic Bohr radius and the screening length due to a finite mediator mass.When this condition for the existence of bound states is met, their binding energy can be written as follows where a 0 = 2/(M α) is the Coulombic Bohr radius (inverse of the Bohr momentum).The discrete energy levels for a Yukawa potential depend on both the principal quantum number n as well as on l.To compute the annihilation width, we project the annihilation terms of pNRY γ 5 onto |n, 0⟩ bound states.In this case, we split the spin singlet and spin triplet bound states.
As one may read from the matching coefficients, nS ≡ |n00⟩ bound states in a spin triplet configuration are absent, whereas for nP ≡ |n1m⟩ states (l = 1) the situation is the opposite.The counterpart of the annihilation cross section eq.(3.24) read as follows for J = 0, 2. R nS and R nP stand for the radial bound-state wave functions.At leading order, their Coulomb limit gives the analytical expressions , whereas in the Yukawa case one has to extract their values at the origin numerically.As argued later in section 3.3, we shall capture bound-state effects by including only the ground state, namely the 1S state.It is therefore useful to define the absolute value of the ground-state binding energy

Bound-state formation and dissociation
Dark matter pairs may annihilate in the form of bound states according to the decay widths given in eq.(3.26).However, bound states have to form in the first place.In order to be quantitative on how such an annihilation process contributes to the overall depletion of dark matter pairs, one has to estimate their formation and dissociation rates.Within the minimal model Lagrangian (2.1) the relevant bound state formation (BSF) process occurs via radiative emission of a single scalar mediator (X X) p → (X X) n + ϕ.The reverse process, namely the bound-state dissociation, is only relevant for the freeze-out in the early universe.Here, a thermal population of scalars can break the bound states if the mediator is sufficiently energetic (roughly of the order of the binding energy of a given bound state).The population of dark scalars is nowadays absent, therefore, bound-state formation is the relevant process for the cross section driving the indirect detection.
In order to compute the bound-state formation cross section, we consider the ultrasoft transitions that are listed in the second line of the pNRY Lagrangian (3.20).More specifically, we are interested in transitions between a scattering state and a bound state, that trigger the formation of a bound state via the radiative emission of a mediator. 6 As outlined in ref. [78], there are various ultrasoft vertices that may contribute to the bound-state formation.The monopole contribution, namely the term at order r 0 in the multiple expansion, vanishes because of the orthogonality of scattering and bound-state wave functions.Then the quadrupole and derivative interactions drive the ultrasoft process (X X) p → (X X) n + ϕ.According to the hierarchy of scales in eq. ( 2.4), we can write the bound-state formation process at finite temperature by factoring out the stimulated emission of the thermal scalar as follows [78] where n B is the Bose-Einstein distribution.When extracting the DM energy density, we will take the thermal average of eq.(3.27) to use it in the effective Boltzmann equation (cf.eq.(4.2)), whereas we take its in-vacuum limit with the appropriate average as discussed in section 5 for present-day annihilations.In the following, the thermal average of a generic cross section is indicated with ⟨σv rel ⟩ and we take Maxwell Boltzmann distributions for the incoming DM fermion and antifermion (for a detailed discussion of the thermal averaging procedure see [90]).The in-vacuum part of the cross section for a given bound state, or exclusive cross section, reads

.28)
6 There can be an additional bound-state formation process in the form of 2 → 2 inelastic scatterings SM1 + (X X)p → SM2 + (X X)n, that are induced by mixing between the scalar ϕ and the Higgs boson (see refs.[53,58,89] for the case of a vector mediator).However, due to the small mixing with the Higgs boson, the typical effective coupling is fairly smaller than α and we neglect this process in this work.Because we assume the trilinear coupling ρ ϕ ϕ 3 = 0 in eq.(2.1), the corresponding process ϕ + (X X)p → ϕ + (X X)n is absent.
where the energy difference between a scattering and bound state is Within pNRY, the BSF cross section factorizes in quantum-mechanical matrix elements, that are often dubbed in the literature as overlapping integrals.In the Coulomb case, it is possible to obtain analytical expressions for such matrix elements [54,56], whereas they have to be numerically evaluated in the Yuakawa case.In this work, we use the wave functions of the scattering states and bound states as extracted with the leading Yukawa potential in eq.(3.22).Because the emitted mediator has a finite mass, the bound-state formation process has a threshold (at vanishing, or very small, relative velocity one has m ϕ < M α 2 γ 2 nl /4n 2 ).We remark that the total bound-state formation cross section entails a sum of all bound states, namely σ BSF v rel = n σ n BSF v rel , where n indicates collectively the quantum numbers of a given bound state (later in section 5 we shall use a different notation for the cross sections, where the quantum numbers n and l will be tracked separately, cf.eq.(5.4)).In our work we shall consider the formation and decays of the ground state only, and estimate the correction of excited states on the relic density, for more details see discussion at the end of section 4.
We specify the general expression in eq.(3.28) for the ground state, and recast the bound-state formation in terms of the parameters ζ and ξ.The phase-space suppression pps nl due to the finite mediator mass is Moreover, we factor out πα 4 /M 2 from the bound-state formation cross section and define a dimensionless quantity that is the analog of the Sommerfeld factors S l ann (ζ, ξ), and we obtain with For illustration, we show S 1S BSF (ζ, ξ) for different values of ξ in the left panel of figure 2, where we fix the scalar coupling to be α = 0.1.The Coulomb limit, which corresponds to ξ → ∞, is also added for comparison.One may see various behaviors for the Yukawa bound-state factors.For ξ = 100 and ξ = 75 in this example, the bound-state factors flatten at large ζ, or equivalently at small velocities, due to the screening of a finite mediator mass.Resonance structures appear also in the bound-state factor S 1S BSF (ζ, ξ) as can be seen in the curve for ξ = 100.For ξ = 20 we notice that the bound-state formation drops severely and BSF ⟩ is performed by taking a Maxwell Boltzmann distribution for the incoming fermion-antifermion pair [39].The bound-state factor is plotted as a function of the time variable z = M/T that is used to solve the Boltzmann equation (4.1).In this case, the Bose enhancement factor [1 + n B (∆E p n )] is included in the thermal average.One can notice how the resonance structure is averaged out and not visible anymore (compare ξ = 100 and ξ = 75 curves).
As for the determination of the DM energy density, one last ingredient is needed, namely the bound-state dissociation.Whenever ionization equilibrium holds, the boundstate dissociation can be simply inferred from the bound-state formation via the principle of detailed balance.For this particular case, one often refers to it as the Milne relation in the literature (see e. g. [52,80] for recent derivations).More specifically the ionization cross section can be written as Here, g X = 2, g ϕ = 1, and g n stand for the d.o.f. of the DM particle, scalar mediator, and a given bound state, respectively.The latter is then inserted in a convolution integral in the momentum of the incoming thermal scalar ϕ to obtain the bound-state dissociation rate where |k| is the momentum of the scalar mediator and Alternatively, one can derive the bound-state dissociation rate from the thermal self-energy of a bound state in pNRY γ 5 [56] (see [91] for the derivation in the context of QCD and heavy quarkonium).For the estimation of the DM energy density, we shall employ the bound-state dissociation of the ground state.
Before moving on, let us comment on the running of the scalar coupling g.The corresponding coupling strength α appears in the hard matching coefficients eq.(3.14)-(3.17), in the Sommerfeld factors through ζ = α/v rel and finally in bound-state formation cross section eq.(3.28) (through the ultrasoft vertex of pNRY γ 5 ).A hard, soft, and ultrasoft scale, respectively, play a role.A running coupling would require tracking the different α(µ) in the various observables.One way to inspect the running of α is to consider the one-loop diagrams (DM fermion and scalar self-energies, and the vertex diagram) and perform the integration by regions [92], i. e. expanding the integrals according to the typical momentum in the loops (soft and ultrasoft) with respect to the mass scales M and m ϕ .Alternatively, one can perform the calculation in the non-relativistic EFTs.Both methods show that no running is induced at scales below M and, therefore, α can be fixed at the hard annihilation scale where it remains frozen. 7

Relic density
In this section, we summarise the main steps and approximations to obtain the DM energy density of the model.This will serve as an input for the interplay between the cosmologically viable parameter space and the corresponding indirect detection prospects.
In this work, we rely on an effective treatment of bound states in terms of semi-classical Boltzmann equations, which allows estimating the bound-state effects on the DM energy density [39,93]. 8In the most general case, the situation is rather complex because there is an equation for the DM particle number density, denoted by n X , and an equation for the number density of each bound state, that results in a set of coupled Boltzmann equations.However, whenever the reactions that govern the rate of change of the bound states are faster than the Hubble rate, the network of Boltzmann equations for the bound states greatly simplifies and turns into algebraic equations [93].In our case, the relevant particle rates are the bound state dissociation rate and the bound state decay width, which are both (much) larger than the Hubble rate for the mass parameters and couplings that we consider in this work.
In doing so, a single Boltzmann equation for n X is found, where the reprocessing of fermion-antifermion pairs into bound states and their decays is accounted for in an effective cross section dn where the Hubble rate can be expressed in terms of the energy density H = 8πe/3/M Pl , M Pl = 1.22 × 10 19 GeV is the Planck mass, e = π 2 T 4 g eff /30, and g eff denotes the effective number of relativistic degrees of freedom.The factor of 1/2 on the right-hand side of eq.(4.1) is needed for DM particles that are not self-conjugated [90].The effective thermally averaged cross section when neglecting bound-to-bound transitions yields where the sum runs over all possible bound states.As already mentioned, we capture the leading bound-state effects by including the ground state only and hence restricting the sum to the 1S state.
In figure 3, left panel, we show the dark matter energy density as a function of the ξ parameter for α 5 /α = 10 −2 , M = 2 TeV and two different choices of the scalar coupling α.The solid curves correspond to ⟨σ eff v rel ⟩ without bound-state effects, namely the second term in eq.(4.2) is dropped.The dashed curves entail both non-perturbative effects.One may see how all curves flatten for large ξ and provide the smallest DM energy density for each choice of the parameters (mass, couplings, and ξ).This is because the Coulomb limit is approached and the non-perturbative effects are maximal.Bound-state effects are absent for ξ ≲ 3 for this choice of the couplings.This is due to the non-existence of the ground state for ξ < 0.84 [76] and to the phase space suppression of the bound-state formation process (cf.eqs.(3.30) and (3.32)), which quickly fades away with increasing ξ values.For the benchmark point α = 0.13 and M = 2 TeV, bound-state effects are crucial to attain a DM energy density compatible with the observed value or, at least, to avoid overclosing the universe.
The right panel of figure 3 offers a complementary perspective.Here, we consider the ratio between the DM energy density as obtained by including or omitting the boundstate formation and decays (i.e. the second term on the right-hand side of eq.(4.2)), in addition to the Sommerfeld factors for the pair annihilation.An interesting trend can be observed.Bound-state formation becomes relatively more important by decreasing α 5 /α ratios, see the orange-dashed, pink-dotted, and dot-dashed-brown curves.This effect is traced back to the different dependence on the scalar and pseudo-scalar couplings of the annihilation cross section of the scattering states σ ann v rel , see eq. (3.24) and the matching coefficients eqs.(3.14)-(3.17),with respect to the bound-state formation cross section σ 1S BSF v rel in eq.(3.31).On the one hand, σ ann v rel depends on α 5 , which triggers important velocity-independent contributions.On the other hand, the bound-state formation is independent of α 5 at the accuracy we have computed such a process in this work (the dominant contribution to the bound-state formation is induced by the scalar coupling for g 5 < g).However, such a trend gets inverted for very small α 5 .This is because the bound-state effects also depend on α 5 through the decay width Γ 1S ann , and the factor Γ 1S ann /(Γ 1S ann + Γ 1S BSD ) entering the effective cross section in eq.(4.2).By taking very small pseudo-scalar couplings, the efficiency factor Γ 1S ann /(Γ 1S ann + Γ 1S BSD ) decreases sensibly and, for α 5 /α = 10 −4 , the bound-state decay is rather suppressed (see the dashed-red curve in figure 3). 9 The relative importance of bound-state effects with respect to the Sommerfeld factors is largely independent of the DM mass (at least in the mass range relevant for this work), that we have fixed to M = 2 TeV in the exemplary case shown in figure 3, right panel.
In our estimation of the pair annihilation via bound states, both for the DM energy density and indirect detection, we include the ground state only, namely 1S= |100⟩.In light of recent findings on the relevance of excited states for DM freeze-out [57,59,94], some comments are in order.These works focus on models with vector mediators from unbroken gauge symmetries.In this class of models, they find a significant logarithmic enhancement of the total bound state formation cross section if a large number of Coulombic states is included [94].This is caused in part by the proliferation of bound states with l ≤ n − 1 for large n.In this work, we are interested in the case of a scalar mediator with a finite mass.As we are dealing with a Yukawa potential, which supports only a finite number of bound states, and not with the Coulombic case, a strong enhancement from summing BSF cross sections up to large n is not expected here 10 .Nevertheless, we have estimated 9 In the limit α5 → 0 the ground state cannot decay and bound-state effects, at the accuracy we are treating them, would be absent, namely Γ 1S ann /(Γ 1S ann + Γ 1S BSD ) → 0 in eq.(4.2).The ground state decay width could be made independent of α5 if one included three-body decays, such as diagrams with three scalar vertices or with a three-scalar-field vertex ρ ϕ ϕ 3 .
10 Requiring the binding energy of the Coulombic bound state to exceed the mediator mass leads to a conservative upper limit on the number bound states that can be produced at all velocities.In the parameter space considered in our phenomenological analysis, this yields n ≤ 22.
the correction from higher states to the effective cross section in eq.(4.2) by including, in the Coulomb limit, the states with n = 2, namely 2S and 2P for four additional bound states (P states with l = 1 have three polarizations of the magnetic quantum number).
Varying the various parameters of the model, i.e. couplings and masses, we found that the corrections to the relic density are always below 10%.The dominant contribution comes from the 2S state because of two features of the model under study.First, the efficiency factor Γ n ann /(Γ n ann + Γ n BSD ) for 2P states is fairly smaller than that of the 2S state.This is due to the different scaling of the corresponding decay width and dissociation rate with the couplings α and α 5 .Second, the ultrasoft vertices of pNRY γ 5 only allow for ∆l = 0 and ∆l = 2 transitions [56,74,78].Thus the ∆l = 1 transitions for 2P states to 2S state induced by the electric-dipole transitions of a vector-like mediator are absent, and 2P states cannot be reprocessed into the more effectively decaying 2S state.A more in-depth study of the tower of bound states for the scalar mediator model goes beyond the scope of our work.

Indirect detection
As discussed previously, the rate of dark matter annihilations in the late universe is too low to affect the total dark matter abundance in a meaningful way.Nevertheless, the effects of such residual annihilations on the flux of cosmic rays can be significant.This opens a powerful method to test annihilating dark matter known as indirect detection.In this work, we are chiefly interested in the ability of current and future gamma-ray telescopes to test our model.Concretely, we will focus on the existing bounds that can be derived from Fermi-LAT observations of dwarf spheroidal galaxies (dSphs) and the prospects for CTA.Note, that due to the limited field of view of Imaging Air Cherenkov Telescopes a different observation strategy is anticipated for CTA, and the best prospects are expected for observations of the center of the Milky Way, Galactic Center (GC) in short.In addition, we also consider bounds on the energy injection around the era of CMB decoupling that have been derived by the Planck collaboration.Thus, we need to model DM annihilations in three distinct astrophysical/cosmological environments: the CMB era, dSphs, and the GC.
All these systems have in common that the typical velocities are much smaller than during freeze-out.However, there still remains a huge variation in their characteristic velocities with the CMB being sensitive to v rel = O(10 −8 ), while dSphs feature v rel = O(10 −5 ) and the Galactic Center even has v rel = O(10 −3 ).An illustration is shown in figure 4, where we depict the DM annihilation and BSF cross sections as a function of velocity for a representative set of parameters.We have highlighted the ranges of velocity that matter for the different environments for comparison.As can be seen, all our cross sections feature strong non-trivial velocity dependencies and the relative importance of annihilations in these systems has to be analyzed with care.Nevertheless, a few qualitative statements can be already made here.First off, since the velocities are significantly smaller than during freeze-out the velocity expansion of the cross sections can be truncated after the leading contribution.Therefore, s-and p-wave refer to the leading piece of the l = 0 and We also shaded the regions of typical DM velocities during freeze-out (FO), at times of the CMB decoupling as well as for dSph and the GC.l = 1 cross sections in the following.Typically, the s-wave cross section for scattering states and the bound-state formation cross section are enhanced in the GC and dSphs compared to freeze-out since the velocities are smaller.Provided m ϕ ≪ v rel M , the same holds for the p-wave process since the Sommerfeld factor compensates for the velocity suppression usually associated with higher partial waves in this regime.Finally, the CMB is chiefly sensitive to the s-wave and bound state cross section since these tend to a constant value for small v rel whereas the p-wave cross section is suppressed at CMB velocities for all masses and couplings considered here.
To make connection with observations one needs to refine the standard treatment to deal with the strong velocity dependence of the cross sections.In the case of the CMB, this is relatively simple since the universe is still highly homogeneous at that point.Thus, the rates can be determined by a simple thermal average as discussed in the previous section.However, the situation is more involved in the GC and in dSphs.Therefore, the following discussion concentrates on the signal at gamma-ray telescopes from these sources.
The key observable in gamma-ray telescopes is the photon flux from the target of interest.The contribution of DM annihilation to the differential photon flux is in general given by (e. g. ref. [95]) dEγ where the sum runs over all possible final states f .The first two integrals run over the observed solid angle Ω and the distance along the line of sight (l.o.s.) ψ.Finally, there are two integrals over the velocities of the two annihilating DM particles.The functions to be integrated consist of the product of the dark matter phase space distribution functions (DF) f X for the two particles and the in general velocity-dependent DM annihilation cross section σv rel .If σv rel is independent of the DM velocity, e. g. for the leading contribution to s-wave annihilation without Sommerfeld enhancement, it can be pulled out of the integrals, and with one obtains the usual J-factor J 0 = dΩ dψ ρ X (r(ψ, Ω)) 2 . (5.3) In the following, we take the dark matter distribution in the dSphs and around the GC to be spherically symmetric and use r = |r| to denote the distance from the center of the object under consideration.When non-perturbative effects are considered (or p-wave annihilations play a role), the velocity dependence of σv rel is more complex and we need to perform the full velocity average.This requires knowledge of the DF.In the following, we assume that we can split up the cross section as σ (l,j) ann v 2(l+j) rel Ionization is negligible at late times such that the effective factor for BSF in eq.(4.2) can be neglected here.This motivates the definition of a generalized J-factor ann (α/v rel , ξ) (5.5) and analogous to BSF.Regarding the annihilation cross section, we are mainly interested in the dominant s-and p-wave contributions which we identify as σ s-wave ann ≡ σ For the rest of the section, we limit ourselves to BSF into the ground state σ BSF ≡ σ (1,0) BSF = πα 4 /M 2 and drop the labels (n, l) in the following.Thus, the DM-induced photon flux can be compactly written as  σ s-wave ann J (0,0) ann,α (ξ) + σ p-wave ann J (1,0) ann,α (ξ) + σ BSF J BSF,α (ξ) .
(5.6) To make further progress here we need to determine the generalized J-factors and the photon number per energy per annihilation.

Distribution function f X
Unfortunately, the distribution function of dark matter in phase space is not directly accessible to observations.Therefore, additional information and/or approximations are required to make progress here.In the following, we use methods for the analysis of galactic dynamics that hold in collisionless systems in a quasi-equilibrium state, see ref. [96] for detailed discussion and derivation, and refs.[95,97,98] for earlier applications to DM halos.An alternative possibility which is not considered here is to sample f X directly from high-resolution numerical simulations of dark matter halos [99,100].
We assume that we are dealing with isolated systems such that the gravitational potential Ψ can be determined from the dark and baryonic matter distribution via a Poisson's equation from the local density where G is Newton's constant and ρ the total matter density.Treating the DM halo as a spherically symmetric system of collisionless particles with density given by ρ X , one can use the Eddington inversion method [96] to obtain a unique and ergodic DF supported by this potential.Here ϵ = Ψ(r) − v 2 /2 is the energy per unit mass of the DM in the gravitational potential.Setting the relative gravitational potential Ψ(r) to 0 at r → ∞ removes the second term.Technically, this boundary condition is unphysical since the gravitational potential of a galaxy is naturally bounded by its neighbors.However, it was shown in ref. [98] that the effect of a large but finite boundary compared to an infinite one is rather small such that we will neglect it.Since Ψ is a monotonic function in r, one can rewrite the Eddington equation directly as a function of r and v [95] (5.9) where r min is given by Ψ(r) − Ψ(r min ) − v 2 2 = 0. Different functional forms for the dark matter halo are considered in the literature.In the following, we will use the NFW-profile [101][102][103] as our ansatz for dSphs.It is given by where r s is the scale radius, which corresponds to d log ρ X d log r r=rs = −2, and ρ 0 fixes the normalization.For our analysis of classical dSphs, we use the best-fit values for r s and ρ 0 as reported in ref. [97].To stay as close as possible to the CTA sensitivity study for the Galactic Center, we model the halo with an Einasto profile in this case [104].It reads with γ = 0.17.The associated gravitational potentials for both profiles can be computed analytically, see e. g. ref. [95].It is more convenient to switch to dimensionless coordinates in the following.We use x ≡ r/r s and the scale mass M s ≡ 4πρ 0 r 3 s h with h ≡ with x min defined analogously to r min and relates to the DF defined in eq. ( 5.9) via To compute the averages over the DFs we normalize them by the DM density and find position-dependent velocity distributions as follows that are normalized to unity, i. e. 4π vesc 0 dv v 2 P r (v) = 4π ṽesc 0 dv ṽ2 P x (ṽ) = 1 where v esc = 2Ψ(r) (and analogously ṽesc ) is the escape velocity.
For the dSphs, baryonic contributions to the gravitational potential can be neglected (Ψ dSph = Ψ X ).However, in the GC baryons contribute significantly to the potential and have to be taken into account (Ψ GC = Ψ X + Ψ b ).We will restrict ourselves in the following to model the baryonic bulge and stellar disk.To apply Eddingtons approach, we need the gravitational potential induced by baryons Ψ b to be spherically symmetric.This is clearly not the case (especially not for the disk).Following ref. [95], we use symmetrized models where M bulge = 1.5 × 10 10 M ⊙ and M disk = 5 × 10 10 M ⊙ are the bulge and disk masses respectively, while b disk = 4 kpc, c 0 = 0.6 kpc are model parameters that encode the spatial extent of these objects.Both are constructed to enclose the same mass as the flattened profiles.

Generalized J s factors
We now turn to the effective J-factors as defined in eq.(5.5).Three of the velocity integrals can be performed directly and only the magnitude of the velocities v 1 and v 2 and the angle ϕ between the annihilating particles remain.Since the integrand only depends on the relative velocity, it is convenient to trade this for the magnitude of the center of mass velocity v cm = |(v 1 + v 2 )/2|, the magnitude of the relative velocity, v rel = |v 1 − v 2 | and z = cos θ, where θ denotes the angle between the two velocities.
Isolating the integration over v rel one gets 2ṽesc 0 dṽ rel P x,rel (ṽ rel )ṽ 2(l+j) rel S (l) ann (Kα/ṽ rel , ξ) (5.17 where we introduced a unit conversion factor K = rsc 2 GMs , and restored factors of c that are implicit in the computation of σv rel in natural units.The integration over ṽcm and z is encapsulated in where the limits of integration are set by the escape velocity and the kinematics with z 0 ≡ ṽ2 esc −ṽ 2 cm −ṽ 2 rel /4 ṽcm ṽrel .In analogy to the thermal average, we define a velocity-averaged Sommerfeld factor.In a slight abuse of notation, we write it as 2ṽesc 0 dṽ rel P x,rel (ṽ rel )ṽ which approaches unity for a decoupled mediator (i.e. no SE) for the s-wave contribution as it should.Thus, we get and analogously for J BSF,α (ξ).This leaves the integration over the spatial coordinates.For dSph, we can safely assume that our solid angle ∆Ω is large enough to cover the entire region where DM annihilation is relevant and exploit that the distance D between us and the source is much bigger than its extent.Thus, we can change coordinates and integrate over a sphere centered around the considered dSph, i. e. dΩ dψ → 4π/D 2 r 2 dr.We consider the classical dwarfs Coma Berenices, Ursa Minor, Draco, and Sergue 1, and parameterize the halos as NFW profiles.The parameters are taken from table 1 of ref. [97].We calculate the enhanced J-factors and present an exemplary set in figure 5 (top).As expected, the J-factors inherit the resonance structure of the Sommerfeld factors and we find a series of peaks once a smooth initial stage has been passed.For large values of ξ the enhancement is very pronounced for all contributions and the hierarchy between the s-wave and the p-wave is reduced considerably.For m ϕ ≲ E b also BSF plays a role and the BSF J-factor rises swiftly with increasing ξ.It clearly exceeds the p-wave factor and is only mildly suppressed compared to the s-wave.Note, however, that due to the different dependence of the split-off prefactors on the couplings and the masses, this information is not sufficient to judge the relative importance of the different contributions for the present-day annihilation rates.
For the GC, we take an Einasto profile and use the halo parameters r s = 20 kpc, ρ s = 0.081 GeV/cm 3 as suggested by a dedicated prospect analysis for CTA [104].We consider the same region of interest (ROI) l, b ∈ [−6 • , 6 • ] where l, b denote the galactic longitude and latitude, respectively, and use D ⊙ = 8.5 kpc as the distance of the Galactic Center to our sun.
We calculate the enhanced J-factors for the Galactic Center and present them in figure 5 (bottom) for the main contributions of s-and p-wave annihilation as well as BSF.The broad structure of the J-factors is similar to the dSph case.However, as expected, the absolute values are larger since the Galactic Center exhibits the largest and closest concentration of the DM.It is also worth noting, that the resonance structure is a bit more washed out due to the higher average velocities compared to the dSphs.With the J-factors taken care of, we can now turn to the next piece needed to determine the gamma-ray flux.

Photon spectrum
The photon spectrum dN γ /dE γ accounts for the number of photons N γ produced by annihilating or decaying DM per photon energy E γ .In our scenario, the annihilations produce pairs of scalars ϕ, which then decay to SM particles via mixing with the Higgs.These SM particles in turn undergo decay and hadronization until only particles that are stable over astrophysical distances remain in the spectrum.Therefore, we need three ingredients to determine the final photon spectrum: i) the branching ratios of ϕ into the various SM final states, ii) the photon spectrum produced by a given SM final state, and iii) a prescription to take into account the large boost of the decaying ϕ's in the galactic rest frame.
In the following, we will collect the relevant branching ratios for the decay of the mediator into pairs of SM particles and determine the photon spectra for cascaded DM annihilation within two mass regimes for the mediator, namely m ϕ ≥ 10 GeV and 2m π 0 ≤ m ϕ ≤ 1 GeV, which we will refer to as the high and low mediator mass range.

Branching ratios for ϕ
The Higgs portal term in eq.(2.2) leads to a mixing between the mediator ϕ and the Higgs.From the point of indirect detection, the absolute value of the mixing is irrelevant, since even very small values lead to a decay of ϕ that is instantaneous on astrophysical scales.Thus, the only question is what ϕ decays into.For m ϕ ≲ 2m h , the branching ratio of the mediator into SM particles equals the one of a Higgs with the same mass.For m ϕ > 2m h the decay channel of the mediator in pairs of Higgs particles opens.The contribution of this channel to the total decay width modifies the branching rations slightly in this regime and we explicitly take it into account.We compute the branching ratios of a hypothetical Higgs with m h ≳ 2 GeV using the expressions compiled in [105].The running of couplings and masses is taken from [106].Above the threshold, the additional decay rate into two : Branching ratios of the mediator ϕ decaying into SM particles as a function of the mediator mass.Assuming the mixing angle between the mediator and the Higgs is small enough, these are equivalent to the Branching ratios of a Higgs with arbitrary mass with the exception of the ϕ → hh channel.
Higgs particles is given by The BR of the ϕ → hh decay becomes comparable to the W and Z contributions above the mass threshold.A display of the branching ratios of all major decay channels within the whole mediator mass range is given in figure 6.For the analysis of the high mediator mass range, we consider the decay channels ϕ → {hh, W + W − , ZZ, gg, tt, bb, cc, τ + τ − }.This method breaks down close to the confinement scale.For m ϕ ≲ 2 GeV we therefore use the branching ratios derived by [107].It is interesting to note here, that in this regime decays into pion pairs play a large role and completely dominate the total width below the kaon threshold.Since π 0 s decay to a photon pair with a branching ratio of 99%, a very large fraction of the energy goes directly to hard gamma rays in this mass range.Hence, we focus on this decay channel here and neglect photons produced in other decays of light mediators.We will only be interested in the continuum flux in setting the limits or prospects and, therefore, this approach is conservative.

Cascade annihilations
We are interested in the photon spectrum produced by DM annihilations.In our case, this process can be broken into two pieces: the annihilation into pairs of ϕ and the decay of ϕ into SM particles.For center of mass energies of the SM pair, i. e. ϕ masses, larger than a few GeV, the second part can be solved with the help of event generators for high-energy physics collision such as Pythia [108] or Herwig [109].Tabulated photon spectra as a function of the center of mass energy E CM for all relevant SM final states have been published in e. g. [110][111][112].In the following we use the results of [111].For m ϕ ≤ a few GeV this approach is unreliable since non-perturbative effects in QCD become large.There have been attempts to derive photon spectra for lighter vector mediators using vector meson dominance combined with data from e + e − to hadrons [113,114].Unfortunately, this is not applicable to scalar mediators and, therefore, we exclude the range 1 GeV ≤ m ϕ ≤ 10 GeV from our analysis.Lighter scalars decay dominantly into mesons and thus the spectra can be constructed from those of meson decays, see below for more details.
This leaves us with the task of translating the photons spectrum computed in the ϕ rest frame to the galactic rest frame.For an isotropically decaying scalar, the spectrum in the center of mass frame of the annihilations is given by [115] where E 0 and E 1 denote the photon energy in the restframe of the decay and the center of mass frame of the annihilation, respectively.The δ function incorporates the Lorentz boost constraint with θ the angle between the direction of the boost and the direction of the photon in the scalar rest frame and Mediators in a mass range 2m π 0 ≤ m ϕ ≤ 1 GeV decay predominantly into pions which makes the derivation of the photon spectrum easy. 11The considered process is X + X → 2ϕ → 4π 0 → 8γ, which is an example of a two-step cascade annihilation and can be dealt with using the method outlined in [117].Since the spectrum of pion decay, dN γ / dx 0 = 2δ(x 0 − 1) with x 0 = 2E 0 /m π 0 , is particularly simple, the cascade spectrum can be computed analytically.In the limit ϵ ϕ,π 0 ≪ 1, where ϵ π 0 = 2m π 0 /m ϕ , the result is where we defined x = E 1 /M .However, these mass ratios are not reached in practice.Therefore, the full result has to be used which reads (5.25) and 0 elsewhere, with (5.26) 11 Even lighter mediators decay into electron-positron pairs and the most stringent limits come from the positron flux measured AMS-02 [16], see e. g. [116] for a DM interpretation.We do not include these final states since we want to focus on gamma rays.
An interesting feature of this spectrum is that it is considerably harder than the power spectra expected for most astrophysical backgrounds.Therefore, the signal is most significant at relatively high energies and has the potential to induce a spectral feature that can boost the sensitivity in a dedicated analysis [118,119].We do not pursue this idea further here and rely on results from a strategy for continuum spectra.

Current bounds and prospects
We now make connections with observations and interpret bounds derived by the Fermi-LAT and Planck collaborations in our scenario.Prospects for CTA are then included in a similar way.
Fermi-LAT limits: We derive our limits from the supplementary material of ref. [10], which is based on six years of Fermi-LAT observations of the classical dwarfs Coma Berenices, Ursa Minor, Draco, and Sergue 1.In addition to the limits for benchmark annihilation channels into pure SM final states, which are not applicable in our case, the collaboration published the energy-bin by energy-bin likelihood as a function of the photon flux [120].This allows to derive approximate likelihood functions for a general DM annihilation spectrum and puts an upper limit on the annihilation rate.We follow the prescription outlined by the Fermi-LAT collaboration and use their data in combination with our photon spectra to derive upper limits on the annihilation rate which can be related to limits on the model parameters using our generalized J-factors.As a check, we also derived limits for the benchmark cases analyzed in [10] following the same procedure.They are in good agreement with the official limits. 12f only one annihilation channel contributes significantly to the photon flux, then we can also place a limit on the individual cross section.An illustration of such a case for a representative set of parameters is shown in figure 7.As can be seen, the limit on the s-wave cross section is rather modest if Sommerfeld enhancement is not included.It lies broadly in the same range as the limits on pure SM final states for these DM masses.The inclusion of the Sommerfeld effect strengthens the bounds by several orders of magnitude everywhere.As expected, the strongest impact is found around the resonances we already observed in the J-factors, which can push the limit even below 10 −29 cm 3 s −1 .A similar analysis for the p-wave case is shown in the lower panel of figure 7. Here, the limits are considerably weaker due to the velocity suppression of the annihilation rate.Without Sommerfeld enhancement the limits are very far from the values relevant for a thermal relic and even with the Sommerfeld effect only the resonantly enhanced parts dip below 10 −23 cm 3 s −1 in our example.The overall limit improves with higher masses since we are keeping m ϕ fixed here which leads to stronger Sommerfeld enhancement as M is increased.
Planck limits: The CMB is sensitive to dark matter annihilations since exotic energy injection can affect the recombination history and change the optical depth of the CMB.
Concretely, the additional energetic particles lead to higher ionization levels at redshifts z ≈ 600 − 1000, which modifies the CMB anisotropies and the polarization [121,122].The most stringent limits on such a non-standard energy injection have been derived by the Planck collaboration [1].Assuming that the redshift dependence of the energy injection rate is controlled entirely by the change of the DM density, which implies that the cross section is velocity independent, they place an upper limit on the combination where f eff is an efficiency factor that encodes how efficiently the released energy is absorbed by the intergalactic medium.We use a conservative estimate of f eff ≃ 0.137 for all channels based on the results of ref. [123].Thus, we find where ⟨σv rel ⟩ is the sum of all thermally averaged cross sections that scale as v 0 rel .In our case, this happens for the s-wave and the BSF cross section since the Sommerfeld and BSF factors have reached a plateau at the relevant velocities (compare figure 4 for an illustration).In contrast, p-wave annihilation is highly suppressed in this regime and does not contribute.The limits on the s-wave cross section for a representative set of parameters are shown in figure 7.As can be seen, the inclusion of the Sommerfeld effect strengthens the Planck limit on the cross section considerably and we find an improvement of several orders of magnitude throughout the mass range.Interestingly, the Planck limit turns out the be relatively similar to the one derived from Fermi-LAT observations.At high DM masses, Planck even outperforms the gamma-ray limits.This can be qualitatively understood based on the scaling of the relevant rates with M .Fermi-LAT is sensitive to the photon flux, which scales as ⟨σv⟩ ρ 2 M 2 , whereas the CMB is sensitive to the energy release rate which scales as M ⟨σv⟩ ρ 2 M 2 and thus has a weaker dependence on the DM mass.
CTA prospects: The dark matter prospects for CTA observations of the Galactic Center have been analyzed carefully by the Cherenkov Telescope Array Consortium in [104].Similar to the Fermi-LAT analysis discussed above, this publication comes with tabulated bin-by-bin likelihoods [124].These can be used to estimate the expected upper limit on the dark matter annihilation rate for any spectrum in a way that is analogous to the one for the Fermi-LAT dwarf limits.An illustration of the prospects that can be derived in this way is presented in figure 7.As expected, CTA has the potential to outperform both Fermi-LAT and Planck limits by orders of magnitude in the range of masses considered here.
Before closing this section, we would like to comment briefly on the dependence of limits and prospects from observations of the galactic center on the halo profile.As the matter in the galactic center is dominated by the baryonic component, the J-factor is not constrained as well as in the case of dSphs and the interpretation of the (prospective) observations is less robust.In particular, cuspy profiled tend to predict larger fluxes compared to cored ones.
In this context, we would like to point out that the Einasto profile does not necessarily fall into the category of a cuspy profile.The inner slope of the profile is controlled by the input parameter γ.With the factor γ = 0.17 adopted in this work, the profile is in-between a cuspy and a cored profile.Nevertheless, the dependence of the prospects on this choice still warrants study.In order to estimate the effect of a more clearly cored profile we again follow the suggestion of the CTA-collaboration and consider an artificially cored Einasto profile [104].Here the density is taken to be constant in the inner part of the halo (below 1 kpc).Note that this ad hoc ansatz does not allow for the application of Eddington inversion and we need to determine the impact of the density profile on the velocity dependent J-factors differently.We compute the usual velocity independent J-factors for s-wave annihilations in both halos and assume that the ratio between the velocity dependent J-factors is the same.Using this procedure we find that the cored J-factor is about 50% smaller which translates to a weaker limit by the same ratio.The comparably small change in the J-factor between different density models in contrast to previous studies has also been addressed in [104] and can be explained by the larger ROI considered here.Naturally, one may wonder how robust this estimate is.In order to check our approach we consider a NFW profile which allows the use of the Eddington inversion method.Computing the ratio of the velocity dependent J-factors, we find that the ratio between the Einasto and the NFW J-factors is largely velocity independent and agrees to 10% with the ratio estimated using the velocity independent (s-wave without Sommerfeld) J-factors.Therefore, we believe that the method outline above is robust.These results are also similar to those of Ref. [125], which studies the CTA sensitivity to Wino and Higgsino dark matter.Comparing an Einasto model to a cored Einasto model (with core size 1 kpc) and using a slightly smaller ROI they find a factor of 3 reduction of the prospect compared to our factor of 2. Finally, we would like to caution, that there are other sources of uncertainties which can affect the limits apart from the J-factors.Taking these into account requires a full analysis of the data and cannot be done based on the likelihood tables provided by [124].Therefore, we do not consider them here, see [104] for a detailed discussion of these effects.

Complementary searches and bounds
In this section we address the connection between the DM model and experimental searches that are complementary to indirect detection.We consider direct detection, BBN limits, thermalization of the dark sector and electric dipole moments.

Direct detection and BBN limits
The mixing between the Higgs and the mediator allows for the elastic scattering of DM on SM particles.The rate of this process is strongly constrained by direct detection experiments (e. g. [28][29][30]) The spin-independent scattering cross section on a nucleon N is given by [72] where m N is the nucleon mass, µ XN the reduced mass of the dark matter-nucleon system and f N ≈ 0.35 the effective coupling of the Higgs to the nucleon 13 .For m ϕ ≪ m h and M ≫ m N this leads to a simple upper limit on the mixing angle Currently the most stringent upper bound on the spin-independent DM nucleus scattering cross section comes from the LZ experiment [29].We always assume that sin δ is small enough not to be in conflict with the direct detection data.As the exact value is not relevant for indirect detection this does not have an effect on the indirect detection phenomenology discussed above.Note, however, that this limit can be combined with limits from BBN to exclude low mass mediators [73].In order to leave the abundances of the primordial elements unaffected one needs to ensure that the mediators decay before the onset of BBN.Therefore we demand that the lifetime of ϕ, τ = 1/Γ ϕ , is shorter than the age of the Universe at the onset of BBN which we take to be T = 1 MeV.This yields As BBN leads to a lower limit on sin δ while direct detection provides an upper one, the combination of this argument with the LZ limits allows to exclude a range of scalar masses.Numerically, we find that an unsuppressed partial width into muons is crucial to ensure a short enough lifetime.Therefore, masses of m ϕ slightly above 2m µ are excluded due to the onset of the kinematic threshold.

Thermal contact between the SM and the dark sector
Our relic density computation assumes that the SM and the dark sector share a common temperature and remain in thermal contact during the freeze-out.While a deviation from this assumption is not a problem as such, it would require a more sophisticated analysis of the freeze-out process and the accuracy of the relic density prediction will vary depending on the details of the thermal decoupling between the sectors.Therefore, it is important to check that our assumptions hold.We follow standard arguments in the literature, see e.g.[127], to estimate the couplings required for efficient thermal contact between the two sectors and we summarize the main steps in the following.
In the symmetric SM phase, the dominant processes are found to be 2 → 2 scatterings of the form ϕH ↔ ϕH and ϕϕ ↔ HH † .As usual, we want the interactions to be efficient and faster than the expansion rate of the universe, which translates to the condition γ 2→2 ≥ 3 H(T ) n ϕ,eq , where γ 2→2 is the rate density (i.e. the right-hand side of the Boltzmann equation) and H(T ) denotes the Hubble rate, which was already defined in sec.4, and n ϕ,eq is the equilibrium number density of the dark scalar.For an order of magnitude estimate, it is sufficient to consider the leading-T behavior of the rates density γ entering the Boltzmann equation, namely γ high−T 2→2 ∝ λ 2 ϕh T 4 /(32π 2 ).By using n ϕ,eq ≈ T 3 /π 2 , one finds the condition to be satisfied for At lower temperatures, the dominant thermal connection between the two sectors is provided by the Higgs-scalar mixing.In this situation the most relevant processes for establishing thermal contact between the SM and the dark sector are f ϕ ↔ f V scattering processes, where f denotes a SM fermion and V is a gauge boson, see e.g.appendix A.2 of [128] for a detailed discussion.The values of sin δ required to achieve thermal contact between the sectors as a function of the freeze-out temperature T fo can be read off from Fig. 2 of the same reference mentioned above [128].In the regime of interest to us, i.e.T fo ≳ 25 GeV, one finds that sin δ ≳ 10 −6 − 10 −7 is required to establish thermal contact.In the bulk of the parameter space considered in this work, these values are comfortably below the experimental limits.Note, however, that for very light scalars with m ϕ ≲ 500 MeV the direct detection constrains become very tight and start to exclude the mixing angles needed to establish thermal contact.This region largely overlaps with the combined direct detection-BBN bound.Figure 9: Two loop diagrams that involve a non-trivial complex phase from the pseudoscalar vertex in the dark sector.Double solid lines stand for the dark matter fermion X, dashed lines the dark scalar ϕ, solid lines the SM fermions, and the red vertex denotes the mixing interaction sin δ f f .

Electric dipole moments
The pseudo-scalar interaction in eq.(2.1) violates CP in the dark sector.Through portal interactions CP violation can be transferred to the SM sector and induce electric dipole moments [129][130][131].In the model of our paper, the mixing with the SM Higgs boson does not produce a pseudo-scalar interaction among the dark scalar ϕ and the SM fermions at leading order, rather one finds only a scalar interaction of the form L int ⊃ − sin δ f f ϕ.Hence, there is no contribution to one-loop topologies.The first contribution to EDM of SM fermions can arise at two-loop level.The relevant diagrams are shown in figure 9, where one scalar and one pseudo-scalar vertex and a DM fermion loop enter.
The most stringent EDM bounds are in place for the electron, namely |d e | < 1.1 × 10 −29 cm e [132].One can estimate the two-loop contribution to the EDM by where we have indicated the possible relevant combinations from the scales running in the loop, the electron mass being the smallest scale in the problem.Even taking the least suppressed contribution from the loop, and by using m e ≃ 0.5 MeV ≃ 2 × 10 10 cm −1 , y e ≃ 2 × 10 −6 , α = 0.1, α 5 = α/10 and sin δ = 10 −2 , we obtain |d e | ≃ 10 −30 cm e, which is smaller than the experimental limit.Note that the value for the mixing angle used in the estimation are larger that those we consider in our study.In addition, the heavy scales in the problem are expected to push the EDM to much smaller values than found in this conservative estimate.Therefore, the electron EDM bound is irrelevant in this model.One may wonder if the situation is more promising if the electric dipole moment of the muon is considered instead but it turns out that the larger Yukawa coupling is not sufficient to compensate for the larger mass and the weaker experimental limit.

Parameter space of thermal dark matter
It is of key interest to understand where the limits and prospects of indirect searches stand compared to the theoretical expectation for thermally produced dark matter.Therefore, we now focus on the parameter space of our model that allows us to explain the full relic density via the freeze-out mechanism.In general, the model is characterized by five parameters: two masses (M and m ϕ ), two interaction strengths (α and α 5 ), and the mixing angle between the mediator ϕ and the SM Higgs.As long as the mixing angle is small enough to avoid limits from direct detection and Higgs decays, its exact value is largely irrelevant to the phenomenology.This reduces the number of effective parameters to four.The relic density is measured with percent level accuracy [1].Therefore, it is possible to reduce the free parameters further by imposing that the correct DM abundance is produced in the early universe.We opt to fix α in this way while varying the two masses and α 5 .With only three parameters left, it is possible to show slices through the parameter space for two of them by fixing the third.As the phenomenology is mainly driven by the masses of the involved particles, we take M and m ϕ as variables in our analysis and show fixed ratios of α 5 /α to explore that direction.
Representative examples of slices through the cosmologically preferred parameter space are shown in figure 10.We have limited the analysis to the region 0.5 TeV ≤ M ≤ 10 TeV and m ϕ ≤ αM , where the upper limit on m ϕ ensures the validity of the NREFT employed in our derivation of the cross sections and induces long-range interactions.Gamma ray constraints and prospects are evaluated for m ϕ ≥ 2m π 0 while the CMB bounds are not restricted by that condition.However, the combined limit by BBN and direct detection searches excludes mediator masses slightly below this threshold (red).Furthermore, we do not consider regions where the relic density computations point towards α ≳ 0.25 (gray).In this regime corrections that are not included in our derivations become large and an analysis of the phenomenology based on our computations can lead to misleading conclusions.These include α and α 5 corrections to the matching coefficients (3.14)-(3.17),and corrections to the binding energy of the bound states.As can be seen, the indirect detection limits are already quite constraining if α 5 is not much smaller than α.This is expected since larger values of α 5 increase the s-wave cross section which is not velocity suppressed at any m ϕ .Together with the Sommerfeld effect, this ensures that essentially all configurations with m ϕ ≲ 3 GeV are excluded by the Planck measurements for α 5 /α = 0.1.At higher mediator masses, CMB exclusions rely on the presence of resonances in the Sommerfeld enhanced J-factors and we find stripes of excluded parameter space that extend up to very high values of M and m ϕ .As the overall cross sections decrease for higher masses these become thinner as the mass increases since the resonance condition needs to be fulfilled more and more precisely.In this regime, the gamma-ray telescopes add more information.As one can see, the structure of the Fermi-LAT limits is comparable to the CMB ones.However, due to the better sensitivity of Fermi, they require a less pronounced Sommerfeld enhancement and, therefore, it has broader stripes that encompass the CMB ones for low M .Towards the high end of the considered mass range, this effect becomes less prominent due to the different scaling of the sensitivity with M as noted previously.Fermi does not add significantly to the CMB bounds for M ≈ 10 TeV.It is very interesting to note, that CTA can be a game changer.With its significantly stronger exclusion prospects, it does not rely on the resonances.This will close the gaps between the stripes and allows to probe almost the whole parameter space.Only very high DM masses and large m ϕ can evade the prospects here.Note that gamma rays also have something to add at low m ϕ .The strip with 2m π 0 ≤ m ϕ ≤ 1 GeV, the low mass-region from now on, is also excluded by Fermi-LAT.These limits are just hidden below the CMB ones to increase the readability of the plot.
For smaller α 5 /α ratios the p-wave and BSF contributions become more important.For α 5 /α = 10 −2 the overall picture remains qualitatively similar; new feature become visible for α 5 /α = 10 −3 .Here the limits from Planck measurements decrease substantially in the low mediator mass regime such that configurations with m ϕ ≲ 1 GeV cannot all be reached by CMB limits alone.However, Fermi-LAT limits are still relatively strong here due to the large number of energetic photons produced in the pion cascade regime.Therefore, most of this parameter space is still ruled out by current experiments.For m ϕ > 10 GeV only quite narrow stripes with a strong resonant enhancement are excluded by present experiments.Prospects for CTA are still excellent and most of the gaps between the enhanced stripes can be closed if its sensitivity fulfills the expectations.Note that the granularity of the gamma-ray limits which becomes noticeable compared to the Planck limits in this figure is a numerical artifact caused by the lower resolution of our Fermi/CTA scans.The determination of the J-factors required for the gamma-ray limits and prospects is computationally intensive and increasing the resolution to the same level as the Planck limits was not feasible with the available computing power.Finally, at α 5 /α = 10 −4 the signals are dominated by p-wave and BSF.Here, existing limits are rather weak and CTA is already needed to close the gaps in the low-mass region.At m ϕ ≥ 10 GeV most of the parameter space is currently unconstrained.CTA is still able to contribute here but for m ϕ /M ≳ 10 −2 , it also relies on resonances in the Sommerfeld factors.Some further strips become visible in this figure.These are driven by the resonances of the pwave Sommerfeld factors which do not coincide with those of the s-wave.All considered, indirect detection is highly effective at testing the model in the parts of parameter space where s-wave annihilation is relevant.In the p-wave and BSF-dominated regime, there are already some interesting limits and good prospects for parts of the mass range.However, a definitive test of the model will remain challenging even with future instruments.

Conclusion
The nature and origin of DM remain two of the most pressing problems in physics today.One particularly attractive possibility is that DM was produced thermally by its interactions with SM particles during the first second after the Big Bang.Despite the simplicity of the basic idea, this scenario allows for a very diverse phenomenology which has inspired a huge range of experimental searches.Among these, indirect searches for the annihilation products of DM in cosmic rays stand out since the relevant observables are directly connected with the processes that drive cosmological production.
Studying such indirect signatures for DM annihilations is particularly timely since new and more sensitive instruments such as the CTA will provide new data in the near future.Therefore, it is important to understand the sensitivity of the CTA to thermal DM and compare its abilities in a fair manner with existing bounds.Since there is a huge range of possible observables, we choose the rather robust limits from Planck and Fermi-LAT observations of dSphs as our reference points.
In this work, we focused on a model of fermionic DM interacting with the SM via a (pseudo)-scalar interaction with a massive scalar mediator.This model is largely unconstrained from direct detection and the LHC, and thus constitutes an interesting benchmark for indirect detection.When the mediator is light compared to the DM mass, special care is required in the derivation of the annihilation cross section.Due to the long-range force mediated by ϕ exchange, non-perturbative corrections commonly known as Sommerfeld enhancement become relevant and the formation of bound states is also possible.Bound fermion-antifermion pairs possess a total energy below that of a free pair and are, therefore, thermally preferred because of a lifting of the Boltzmann suppression.Hence, they can have an important effect on the relic density even though the formation cross section is suppressed by additional powers of the coupling.Both Sommerfeld enhancement and BSF are of crucial importance for light mediators and can drastically alter the relevant rates in the early universe and astrophysical environments.We employ NREFT and pNREFT techniques, which allow us to take these effects into account systematically in the computation of the relevant cross section, and we derive the DM relic density with the usual Boltzmann equation method.
As a consequence of the non-perturbative effects, the cross sections of interest exhibit a strong and non-trivial dependence on the relative DM velocity.This breaks the common lore of indirect detection that only the leading velocity-independent contribution to s-wave annihilation is detectable in a realistic experiment.Therefore, we employ an extended version of the standard formalism to compute the appropriate velocity averaged J-factors.This requires position-dependent velocity distributions for the targets under consideration.We derive these using Eddington inversion from the density profiles for four dSphs and the Galactic Center.This allows us to predict the photon flux from DM annihilation from the particle physics parameters and impose upper limits (identify the prospects).
Combining this with our predictions for the relic density allows us to assess the status of thermal DM in this model.If the suppression of α 5 relative to α is modest, s-wave annihilations play a large role and there are strong constraints on light mediators irrespective of the DM mass.At higher m ϕ the resonant enhancement due to the Yukawa Sommerfeld factors becomes important and allows to exclude the regions close to the resonance with Planck and Fermi-LAT observations.Interestingly, the gaps that remain in the testable parameter space can be closed by CTA in this case.Testing the model becomes more difficult for smaller α 5 /α.Here, the sensitivity of Planck and Fermi to the cosmologically preferred parameter space diminishes and CTA is crucial in this regime since it will allow testing large regions of the parameter space for the first time.

γ 5 pNRY γ 5 √Figure 1 :
Figure 1: Hierarchy of energy scales and effective field theories considered in this work for the DM model Lagrangian density defined in eq.(2.1).In-vacuum and thermal scales are shown.

Figure 2 :
Figure 2: (Left) Bound-state formation factor of eq.(3.32) as a function of ζ for ξ = 20, 75, 100.The Coulomb limit is shown for comparison.(Right) Thermal average of the bound-state formation factor with the inclusion of the Bose enhancement [1 + n B (∆E p n )] for the same choice of the ξ parameter as in the left panel.

Figure 4 :
Figure4: Leading order contributions to the DM annihilation cross section from s-wave and p-wave annihilation as well as BSF as a function of v rel for fixed M = 1 TeV, m ϕ = 1 GeV, α = 0.1 and different pseudo-scalar couplings α 5 (note, that for the p-wave / BSF cross section, a change in α 5 plays a subdominant / no role and is therefore not depicted).We also shaded the regions of typical DM velocities during freeze-out (FO), at times of the CMB decoupling as well as for dSph and the GC.
) for Dirac DM.Let us go through this expression piece by piece.The object dNγ dEγ describes the number of photons per energy per annihilation.If more annihilation channels are open it can be constructed from the photon numbers per energy of the different channels weighted by the branching fractions B f into the channels, i. e. dNγ dEγ = f B f dN (f ) γ

Figure 8 :
Figure8: Contours for the condition γ 2→2 /(3Hn ϕ,eq ) = 1.The dashed-orange and solidmagenta curves stand for the condition thermally averaged cross section as with full account of the scalar mass m ϕ (smallest and largest values respectively); the dotted-brown line stands for the estimation based on dimensional analysis only.
)where g eff (T ) denotes the effective number of SM degrees of freedom.For T ∈ [150,10 3 ] GeV one finds λ ϕh ≃ [10 −7 , 5 × 10 −7 ].A more accurate extraction based on full thermally averaged cross sections that accounts for the finite mass of the dark scalar leads to the results shown in Fig.8.The more efficient elastic scattering ϕH ↔ ϕH allows for λ ϕh ≃ 10 −6 for the whole range of scalar masses considered in our study.

Figure 10 :
Figure 10: Slices through the cosmologically preferred parameter space in the M -m ϕ /M plane for four different representative values of α 5 /α (10 −1 , 10 −2 , 10 −3 and 10 −4 ).The value of α is fixed by the relic density requirement.The regions excluded by Planck are shown in turquoise and the Fermi-LAT exclusion in gold.Prospects for CTA are shown in pink.The red region depicts mediator masses, which are excluded by the combined BBN and direct detection limit on the mixing angle.In the gray region the relic density requires α > 0.25.The black lines indicate the characteristic values of m ϕ that delineate the different gamma ray production regimes.