Cosmic-ray antiproton constraints on light singlino-like dark matter candidates

The CoGeNT experiment, dedicated to direct detection of dark matter, has recently released excess events that could be interpreted as elastic collisions of $\sim$10 GeV dark matter particles, which might simultaneously explain the still mysterious DAMA/LIBRA modulation signals, while in conflict with results from other experiments such as CDMS, XENON-100 and SIMPLE. It was shown that 5-15 GeV singlino-like dark matter candidates arising in singlet extensions of minimal supersymmetric scenarios can fit these data; annihilation then mostly proceeds into light singlet-dominated Higgs (pseudo)scalar fields. We develop an effective Lagrangian approach to confront these models with the existing data on cosmic-ray antiprotons, including the latest PAMELA data. Focusing on a parameter space consistent with the CoGeNT region, we show that the predicted antiproton flux is generically in tension with the data whenever the produced (pseudo)scalars can decay into quarks energetic enough to produce antiprotons, provided the annihilation S-wave is significant at freeze out in the early universe. In this regime, a bound on the singlino annihilation cross section is obtained, $\sigv\lesssim 10^{-26}\,{\rm cm^3/s}$, assuming a dynamically constrained halo density profile with a local value of $\rho_\odot = 0.4\,{\rm GeV/cm^3}$. Finally, we provide indications on how PAMELA or AMS-02 could further constrain or detect those configurations producing antiprotons which are not yet excluded.


Introduction
The CoGeNT [1] Collaboration has recently claimed for excess events in their data, while not confirmed by complementary searches like XENON-100 [2,3] or SIMPLE [4]. Although CDMS-II had reported a few excess events [5], while not statistically significant, their latest results are also in conflict with the CoGeNT data [6]. The situation has become even more confusing, though quite exciting, since CoGeNT announced hints for annual modulation of the event counting rate [7]. Scatterings of weakly interacting massive particles (WIMPs) with masses around 10 GeV off the detector targets could explain these measurements, and interpretations in some specific particle physics models have been discussed in e.g. [8] and [9] (see also [10] for less conventional models). Intriguingly, this mass range turns out to be close to what is needed to fit the still mysterious DAMA/LIBRA signals [11,12,13]. Nevertheless, the authors of Ref. [9] have notably shown that some of the parameter space individually favored by the cited experiments may not be compatible among each other (see also [14,15,16]).
Email addresses: davidg.cerdeno@uam.es (David G. Cerdeño), timur.delahaye@uam.es (Timur Delahaye), lavalle@in2p3.fr (Julien Lavalle) 1 Multidark fellow These controversial data have triggered a lot of activity in the field, and different phenomenological directions have been pursued. While light neutralinos arising from the minimal supersymmetric extension of the standard model (MSSM) are still debated (see arguments against this proposal in e.g. [17], to which very detailed highlights are opposed in [18]), singlet extensions of supersymmetric scenarios (including the nextto-MSSM -NMSSM) have received a particular attention, e.g. [19,20,21,22,23]. The NMSSM, which provides a very elegant way to solve the so-called µ-problem occurring in the MSSM (see [24,25] for reviews), could actually lead to a phenomenology consistent with the direct detection hints, though not generically, and some interesting regions of the parameter space are worth being investigated more deeply [26,21,22]. In slight contrast, more general singlet extensions uncorrelated with the µ-problem can quite easily be tuned to accommodate the data, e.g. [23], because more flexible. In these models, the dark matter candidate is usually a light singlino-dominated neutralino.
Interesting constraints on light WIMPs may actually come from colliders [27,28,29], but most of the singlet models mentioned above seem to escape them [21,22]. From the astrophysical point of view, the annihilation of such light WIMPs may also generate gamma-ray fluxes at the edge of being excluded with current measurements [30]; observational limits on highenergy neutrinos from the Sun also give interesting constraints, Preprint submitted to Nuclear Physics B February 1, 2013 as shown in e.g. [31]. Another very powerful astrophysical constraint comes from measurements of local cosmic-ray antiprotons, which severely limits direct annihilation into quarks [32].
In this paper, we aim at extending the latter analysis to singletlike phenomenologies in which singlino-like dark matter particles χ annihilate into light scalar and/or pseudo-scalar particles (generically referred to h and a in the following), which can further decay, for instance into quarks. For the sake of generality, we will adopt a model-independent approach by using an effective Lagrangian relevant to any particular singlet model, and will derive constraints on all involved masses and couplings. We note that indirect detection constraints come into play only when annihilation is not velocity-suppressed, which implies that the dominant annihilation channel is χ χ → a h. This defines the mass range we are going to focus on, namely 2 m χ ≥ m a + m h . The outline of the paper is the following. In Sect. 2, we present our effective Lagrangian approach and discuss the annihilation cross section and issues related to the relic density. In Sect. 3, we briefly review the calculation of the spinindependent WIMP-nucleus cross section, and delineate the parameter space to be considered for interpretations of the direct detection hints. In Sect. 4, we derive the cosmic-ray antiproton spectra generated by decays of light scalar and pseudo-scalar Higgs bosons, and we determine the constraints current flux measurements put on the effective parameter space. We conclude in Sect. 5. Further details on our calculations are given in the appendices.

Effective Lagrangian approach to singlino-like dark matter scenarios
The phenomenology of light singlino-like dark matter scenarios is mostly set by (i) the couplings between the neutralino, a Majorana fermion, and the light scalar and pseudo-scalar particles of the Higgs sector, (ii) the masses of the involved fields, and (iii) the couplings of the light scalars to standard model particles. We address the former point here; the latter will be discussed in Sect. 2.2. Denoting χ the neutralino field and φ i a generic field that contains a scalar part h i and a pseudo-scalar part a i such that φ i ≡ h i + i a i , we consider the following effective interaction Lagrangian: where we λ i, j,k is a dimensional real coupling (η i jk is a combinatory factor), and where C χi ≡ c χi + ic χi γ 5 =C χi (2) features both the couplings of the scalar and pseudo-scalar fields φ i to the dark matter fermionic field χ. We further define c χi (scalar coupling) andc χi (pseudo-scalar coupling) as real numbers.

Higgs sector
We consider singlet extensions of the MSSM in which the lightest neutralino is singlino dominated, and the light Higgs scalar and pseudo-scalar are singlet dominated. Nevertheless, the physical Higgs fields are still mixtures of the Higgs doublets H u and H d and of the singlet field S , which allows us to define the mixing matrix S as where s k and c k stand for sin k and cos k, respectively. The mass eigenstates corresponding to the lightest pseudo-scalar and scalar fields, which we respectively denote a and h, are conventionally given by the first row: where A MSSM corresponds to the same linear combination of the imaginary parts of H u and H d that describes the pseudo-scalar Higgs in the MSSM. The mixing angles can be determined in specific models from the mass terms and couplings featuring the superpotential and the soft potential. The Higgs mixing will actually play an important role not only when calculating the direct detection signals, as discussed in Sect. 3, but also to determine the decay channels of the light Higgs fields, which the antiproton production depends on (see Sect. 4.1). It also helps us understanding the collider constraints, as discussed later on.
In this paper, we are assuming the lightest CP-even Higgs to be singlet-like, in order to avoid collider constraints, and therefore we take c ϕ ≥ 0.9. For the two heavy CP-even Higgs bosons, we adopt a MSSM-like scenario in which the lightest one is of u-type (H u ) and the heavy one is approximately H d . Thus, we fix c θ = c ψ = 0.9.
Last but not least, as it is customary in the MSSM, we define tan β as the ratio of the Higgs doublet fields vacuum expectation values, tan β ≡ H u / H d .

Couplings to standard model fermions
To complement this effective setup, we also need to define the couplings between the scalar and pseudo-scalar fields to the standard model (SM) fermions, which will further determine the branching ratios of their decay products (annihilation final states) and the strength of the WIMP-nucleon interaction (direct detection signals). For a field φ i having both a scalar and a pseudo-scalar part, the couplings to SM fermion fields differs for up-type quarks (or neutrinos) and down-type quarks (or charged leptons) as: where m u (m d ) is the up-type (down-type) quark mass.
These couplings are useful to compute the decay branching ratios of the light Higgs fields, which will further determine the induced antiproton spectrum. In practice, we use the light (pseudo-)scalar decay widths that we derive as: where N c is a color factor (3 for quarks, 1 for leptons), and the latter expression is relevant to the cases where the process h → a a adds up to decays into SM fermions (m h ≥ 2 m a ).
In the following, we will only consider those interactions between singlinos and SM particles which are mediated by the scalar and pseudo-scalar Higgs bosons, assuming all other contributions to be negligible. In other words, we assume that other sectors of the complete underlying theory decouple. This is a reasonable assumption given the current mass limits on squarks [33,34] and charged sleptons [35]. Moreover, as regards the Higgs sector, we impose that the interactions with standard matter are driven by light field exchanges only (see Sect. 3).

Some collider constraints
It is important to derive ranges for our effective model parameters, namely masses, mixing angles and couplings, which are consistent with current collider constraints. These constraints necessarily involve the couplings to SM particles, which have been made explicit in Eq. (5). We basically use the constraints provided in [23] for light scalars, and in [36] for light pseudoscalars. The only additional feature we impose to our singlino dark matter candidate is to reproduce the region compatible with the WIMP interpretation of the CoGeNT signal, given the parameter space delineated for the scalar and pseudo-scalar fields below.
Limits can be derived because of the potential contributions light scalars and pseudo-scalars may give to rare mesons decays. From measurements of B decays by Belle [37] and SM background calculations by [38], the following limit has been derived by the authors of [23] for the light scalar mass, given its non-singlet content s ϕ : 0.058 The decay branching ratio Br(h → µ + +µ − ) is readily computed from Eq. (5). This roughly sets a lower limit for m h in terms of the bottom quark mass, namely m h m b . There is an additional and independent limit on s ϕ which can be derived from LEP data on Higgs searches, especially from the process e + e − → Z h → νν h [39]. This yields the upper limit s ϕ 0.08, valid as long as m h < 2 m a , i.e. as long as the decay channel h → a a is forbidden [23].
Constraints on the singlet content of the pseudo-scalar can be extracted from e.g. the CLEO Higgs search data from the Υ(1S ) decays [40]. They are expressed in terms of the reduced coupling X d of the light pseudo-scalar a to down-type quarks and leptons, which is normalized with respect to the coupling of the CP-even SM Higgs boson, and a very conservative limit is given by: as read off from Fig. 1 of Ref. [36]. This is definitely not a stringent limit for the phenomenology under scrutiny here, since this has no consequence on the pseudo-scalar decay products (the mixing angles factorize out in the decay branching ratios). This still ensures that the singlino annihilation into SM fermions via light pseudo-scalar exchange is always suppressed with respect to the annihilation into light (pseudo-)scalars. Finally, we set a lower limit of m a ≥ 1 GeV in the following to avoid relying on too large fine-tuning zones in more theoretically constrained parameter spaces.

Relic density
One of the most important constraints on dark matter models obviously comes from requiring that they lead to the correct cosmological abundance, which usually severely reduces the available parameter space. Despite the uncertainties affecting our knowledge of the expansion rate of the universe before the big-bang nucleosynthesis (BBN) (see e.g. [41]), the general formalism to derive the relic abundance, based on solving a covariant Boltzmann equation, has been established and tested for a long time [42,43,44,45,46,47]. There are some public codes like DarkSUSY [48] or MicrOMEGAs [49], which numerically solve the Boltzmann equation, originally designed for supersymmetric models. Both these codes essentially make use of the same numerical scheme and of the same tables for the effective degrees of freedom, which rely on earlier works [44,46,47]. Consequently, differences between them are expected to come from the annihilation cross section calculation only. An alternative to these numerical tools, but at the cost of loose accuracy, is to use approximate analytical expressions [45,46,50]. Here, we use our own numerical tools, more detailed in Appendix B to which we refer the reader for the definitions of the relevant parameters. Denoting Y > f the comoving dark matter density after decoupling such that Y > f Y eq , its value today Y 0 can be calculated from the following expression: where x = m χ /T (T is the thermal plasma temperature). Note that the second term in the right-hand side, featuring the annihilation cross section σv (see Sect. 2.5) and the effective degrees of freedom g , is usually dominant. A numerical code is basically useful to compute σv and the effective degrees of freedom g . The relic abundance then merely reads: where we have translated the critical density ρ c in terms of the Hubble parameter H, using the relation H 2 = (100 × h) 2 = 8 π ρ c /3. The parameter s 0 2891.23 cm −3 is the entropy density today -we have taken a CMB temperature of 2.72548 • K [51]. One of the key issues in calculating the relic density of ∼ 10 GeV WIMPs is the QCD phase transition. Indeed, these WIMPs typically decouple at temperatures T f ∼ m χ /20 ∼ 400 MeV, when quarks (and gluons) go from asymptotic freedom to confinement into hadrons. This transition is therefore characterized by an important change in the relativistic degrees of freedom g appearing in Eq. (9), which can be a source of systematic errors in the predictions (as was already clear from e.g. Figs. 5 and 6 of [44]). For instance, slightly modifying the phase transition temperature T QCD can lead to significant changes in the final relic abundance, by up to 50% (a factor of ∼ 4 in g ), for those models which happen to decouple after instead of before the transition, after modification of T QCD (see Eq. 9). The point is that we do not know much about this phase transition yet despite the recent progresses obtained in lattice QCD [52], and theoretical uncertainties in this regime might still be sizable [53]. The current state of the art [54] is implemented in e.g. DarkSUSY, but modifications of the transition temperature is not user-free.
To try to bracket potential theoretical uncertainties coming from this phase transition, we have recomputed the effective degrees of freedom g for different values of T QCD , considering a close to first order phase transition. An example is shown in the left panel Fig. 1, where we display g eff 1/2 , h eff 1/2 and g 1/2 for two different assumptions T QCD = 150 or 400 MeV (see Eq. B.2 in Appendix B for the definitions). The right panel exhibits the consequence on the dark matter relic abundance, assuming a constant annihilation cross section set to its canonical value σv = 3 × 10 −26 cm 3 s −1 , and only varying the WIMP mass m χ . We clearly see that in the region of interest, m χ ∼ 10 GeV, the relic density can decrease or increase by a factor of ∼ 2, depending on whether decoupling occurs just before or just after the QCD phase transition -this can be understood from Eq. (9), where we see that Y 0 ∼ ∝ (g 1/2 σv ) −1 . Therefore, some models in this mass range can easily be cosmologically excluded for one choice of T QCD but allowed for other choices. In the plot, as will be the case in the rest of the paper, we have used the WMAP-7 data for the present dark matter abundance, i.e. Ω χ h 2 = 0.1109 ± 0.0056 [55], which appears as the green band.
This illustrates why uncertainties on the quark-hadron transition and on T QCD should be taken into account, as often underlined in the past, like in [56]. A detailed implementation of this transition is beyond the scope of this paper. Still, we consider the two extreme cases depicted in the left panel of Fig. 1, which is likely sufficient to encompass the full theoretical error. Finally, we stress that the impact on indirect detection signals amounts to the same factor of ∼ 2. The annihilation cross section needs to be relatively larger to get the correct abundance when dark matter particles decouple after the QCD phase transition, as clearly seen in the right panel of Fig. 1. We have checked that using other predictions for the evolution of the relativistic degrees of freedom and their change during the QCD phase transition, for instance those included in the DarkSUSY package, led to a 10% difference in the relic density predictions, or equivalently in the annihilation cross section 2 . Therefore, taking a 3-σ band for Ω χ h 2 is well enough to encompass the theoretical uncertainties in the effective degrees of freedom for a given T QCD . Thus, we will consider any configuration as cosmologically relevant whenever leading to a relic density of Ω χ h 2 = 0.1109 ± 0.0168.

Annihilation cross section
The annihilation cross section is relevant not only to indirect detection, but also to direct detection, since in both cases it determines the cosmological abundance today. In the following, we assume that the only particles involved in the annihilation process are the light scalar and pseudo-scalar Higgs bosons, h and a respectively. This is justified by the fact that we are in the very light mass and singlet/singlino dominated regimeother Higgs bosons, charginos and sfermions have much larger masses, while similar couplings. In such a regime, annihilation implies the following final states: χ χ → a a, χ χ → h h and χ χ → a h, which can be summarized thanks to the Feynman diagrams represented in Fig. 2.
As it is well known (e.g. [57]), χ χ → a h is the only treelevel process among the three above possibilities to have a nonvanishing zero-velocity limit, i.e. the only process relevant to indirect detection. It is kinematically allowed as long as 2 m χ ≥ m a + m h , which defines the mass range we are dealing with in this analysis. The full calculation of the annihilation cross section is detailed in Appendix A. Here, we just summarize a few helpful results. It usually proves interesting to expand the thermally averaged annihilation cross section as a power series of the inverse of the variable x ≡ m χ /T ∝ v −2 (x f ∼ 20 at freeze out, and x 0 ∼ 10 7 in the Galaxy today): Both coefficients a and b take analytical expressions in the present study. The S-wave contribution is given by (we neglect the widths of the exchanged particles 3 ): 2 We have also made some comparisons between our calculation and the one provided in MicrOMEGAs, and found an agreement within ∼ 1% when using the MicrOMEGAS table for the degrees of freedom (identical to the DarkSUSY  table). 3 Particle widths are included in the full numerical calculation.
This equation fully sets the annihilation cross section for the indirect detection signals that we will consider later. We immediately see that there is a competition between the s and t/u diagrams through the couplings λ aah and c χχh . Note that the latter is also the one that fixes the contributions of CP-even Higgs bosons to the direct detection signal (more details in Sect. 3), which sketches a way to link the indirect detection constraints to the direct detection signal: the relation is straightforward in the case λ aah = 0, corresponding not only to the suppression of the s-channel of χ χ → a h but also of that of process χ χ → a a. Finally, for λ aah ∼ c χχh and close-to-degenerate masses, the t/u diagram dominates the cross section, which again induces a significant correlation between the direct detection signal and the indirect detection constraints.
We also remark that the coupling c χχa does drive the a-term amplitude, which is actually quite obvious since χχ is bound to be a pseudo-scalar state at null velocity. Suppressing this coupling is a very simple way to avoid indirect detection constraints. In that case, annihilation before freeze-out would proceed mostly through χ χ → h h, but also partly through χ χ → a a via the s-channel, provided these final states are kinematically allowed.
Although the connection between the elastic interaction with nuclei and the annihilation into a h may be understood from the above equation, we still have to make sure that the a-term of the cross section expansion is the dominant one in the early universe, in order for it to be fixed by the relic abundance. Indeed, should the processes χχ → aa, hh be dominant at freeze out, then that would lessen the bounding strength of indirect detection and also make us lose potential correlations between direct and indirect signals. We recall that the correlation is definitely lost when m χ ≥ m h (m a ) and 2m χ < m a + m h , for which χχ → hh (aa) is the only annihilation possibility left; we do not consider such cases here.
For completeness, we provide the first order terms b aa and b hh associated with annihilation into a a and h h, which drive the relic density when 2 m χ < m a + m h (still neglecting the particle widths): 16 These terms are the leading terms of the annihilation cross section when 2 m χ < m a + m h . We emphasize that even in the case 2 m χ ≥ m a + m h , i.e. when a ah is not suppressed, the velocitydependent terms can not only partly contribute to the full annihilation cross section at decoupling in the early universe (up to ∼10-20% when couplings and masses are taken similar, and considering also b ah , not given here), but can also dominate in some cases, e.g. when the couplingc χχa is suppressed. This means that the full annihilation cross section needs to be considered for relic density calculations, which we have done in this study. In fact, we did not use this Taylor expansion approach to compute the thermal average of σv in the relic density calculation, but the exact numerical method presented in [46,47]. Still, the velocity independent term a ah given in Eq. (12) fully characterizes the annihilation cross section relevant to indirect detection. Moreover, Eqs. (13) and (14) will prove very useful to interpret the antiproton constraints (see Sect. 4).

Direct detection of singlino-like particles
In this paper, we aim at confronting the light singlino parameter space compatible with the CoGeNT region to the antiproton constraints. It is therefore important to provide some details on how the CoGeNT region can be covered in the frame of our effective approach.
The singlino couplings to the Higgs mass eigenstates will be denoted as c χ j , in consistency with the Lagrangian introduced in Eq. (1).
Elastic collision of a singlino off a target nucleus would take place through the exchange of a CP-even Higgs along a tchannel, as depicted in Fig. 3. This scattering can be described in terms of an effective Lagrangian [57], where the coupling α q i results from adding the contributions of the exchanges of the three CP-even Higgs bosons and reads, where i = 1, 2 for up(down)-type quarks, B 1 = sin β and B 2 = cos β. The mixing matrix S has been defined in Eq. (3). Note that the term of the above sum corresponding to our light Higgs boson h can be straightforwardly rewritten as α h q i = c χ j C hqq /m 2 h , where the coupling C hqq has been defined in Eq. (5).
The total spin-independent WIMP-proton scattering cross section yields where m p is the proton mass and The quantities f p T q i and f p TG are the hadronic matrix elements which parameterize the quark content of the proton (the second term is due to the interaction of the WIMP with the gluon scalar density in the nucleon). We will use f p u = 0.023, f p d = 0.033 and f p s = 0.26 which follows from the latest results on the pionnucleon sigma-term [58], which affects the determination of the strange content of the quark.
The leading term in the scattering cross section is usually due to the coupling to the s quark (the contribution from quarks u and d being almost negligible and that from f p TG being at least a factor 5 smaller in the determination of f p ). Under this assumption it is possible to derive from Eq. (16) a condition for the lightest Higgs contribution to be the dominant one in terms of the Higgs masses and the various couplings. In particular, we will demand c χ1 and c χ1 Notice that since we are considering very small masses for the lightest Higgs boson, these conditions are normally easy to fulfil. In our case they will imply a constraint on the Higgs doublet contribution and on the neutralino-Higgs coupling. For example, for a typical configuration with m H o 1 = 10 GeV and m H o 2 = 120 GeV the above condition implies c χ1 S 1 1 > 7 × 10 −2 c χ2 S 1 2 . Under these conditions the spin-independent cross-section can be approximately expressed as In order to explore the region consistent with the latest hints for very light WIMPs in direct detection experiments [12,1,7], we have constrained the input parameters to produce a singlinoproton spin-independent cross section between 3 × 10 −5 pb and 3 × 10 −4 pb, and singlino masses in the range 6 − 15 GeV. This window in the (σ S I χ−p , m χ ) plane contains the region compatible with the first observation by the CoGeNT collaboration [1] and is extended in order to account for variations of this region when astrophysical uncertainties are taken into account [9,59] and to include the area consistent with DAMA/LIBRA result if interpreted in terms of WIMP recoils on NaI (assuming negligible channeling). Of course, the region consistent with the latest data release by CoGeNT [7] is also contained in the whole area.
To try to delineate a parameter space as complete as possible, we have performed different random scans over the full parameter space defined by the involved free masses and couplings, taking different values of tan β for each. Since we require a given range of spin-independent cross-section and, as we said above, we consider only those cases where the lightest Higgs contribution to σ SI χ−p is dominant, this implies a constraint on the lightest Higgs composition and its coupling to singlinos. In particular, for m χ = 6 − 15 GeV, Eq. (21) leads to We note that this condition implies that for large values of tan β the necessary value of the coupling c χ1 in order to reproduce the CoGeNT data can be smaller. In practice, the various input parameters have been randomly drawn in the following ranges: where the scan in the couplings is performed in the logarithmic scale. Note that the lower bound on the couplings λ is not a fixed number, but instead twice the coupling of any light boson to b quarks (see Eq. 5). This is to ensure that (the s-channel) annihilation dominantly proceeds into 2-body combinations of bosons a and h, the pair production of SM fermions being then suppressed and irrelevant in the computation of the relic density as well as of the indirect detection signals. This will induce a correlation between couplings λ and C χh through the range imposed on the spin-independent cross section [see Eqs. (5), (16), and (21)]. The lightest CP-even and CP-odd Higgs masses are assigned a random value, imposing the constraint m h + m a ≤ 2 m χ in order to guarantee that annihilation into a pair h, a is possible. We have imposed the experimental constraints on the Higgs masses and couplings as detailed in Sect. 2.3 and the condition that the singlino-proton spin-independent cross section falls within the window σ SI We have thus kept those configurations which fall in a region that contains the area compatible with the CoGeNT excess. As regards tan β, after having tried some random scans, we have selected two extreme values only, which are sufficient to fully illustrate the impact of this parameter. In Fig. 4, we show our scan results in terms of σ SI as a function of m χ for both values of tan β, and for two different QCD phase transition, with T QCD = 150 and 400 MeV. The color index refers to relic density and antiproton constraints, and will be specified and discussed in Sect. 4.3. All samples we consider in the following are made of ∼ 10 5 random models, such as those reported in Fig. 4.
To conclude this section, we must emphasize that some uncertainties in the hadronic matrix elements may affect the predictions for the direct detection rate of singlinos. The largest effect is due to the indetermination in the strange quark content of nucleons. This propagates into the theoretical predictions for the spin-independent cross section and can be responsible for a variation of about an order of magnitude [60,61,62,63]. This induces an uncertainty in the values of the couplings for which the CoGeNT region can be reproduced and therefore enhances the range of these with respect to those plotted in Fig. 4. A recent illustration has been given in e.g. Fig. 5 of Ref. [64]. We will further discuss this issue in Sect. 4.3.2.

Constraints from local cosmic-ray antiprotons
Cosmic-ray antiprotons were proposed as a potential tracers of dark matter annihilation in [65]. Dark matter-induced antiprotons have been studied in the context of supersymmetry [66,67,68,69,70,71] and extra-dimensions [72,71], and related signals are likely very difficult to observe. Current observations by e.g. the PAMELA [73] or BESS [74,75,76,77] experiments are in agreement with the predictions of secondary antiprotons [78,79,73], i.e. those created from nuclear interactions of standard astrophysical cosmic rays with the interstellar gas, although some models still fail to saturate the data, leaving a small room for exotic primary contributions [80,81]. The case for light supersymmetric neutralinos was studied in [82].
Based on recent studies on cosmic-ray transport parameters [83] and on the local dark matter density [84,85], it was shown that ∼ 10 GeV light WIMPs annihilating into quark pairs were in serious tension with the low energy cosmic-ray antiproton data [32]. In this paper, we are focusing on a quite different phenomenology where WIMPs annihilate into light scalar h and pseudo-scalar a Higgs bosons, which may further decay into quarks. These light bosons do not have necessarily the same masses, so one can for instance be more Lorentz boosted than the other. Therefore, the energy distribution of the potentially created antiprotons is expected to differ from those arising from direct annihilation into quark pairs, and the results obtained in [32] can hardly be extrapolated. In this section, we first determine the antiproton energy spectrum associated with the decay of light Higgs bosons into quark pairs before calculating the dark matter-induced flux at the Earth.

Antiproton production
Before estimating the antiproton flux at the Earth, one needs to know the antiproton spectrum before propagation, i.e. the one arising from the decay/hadronization of the dark matter annihilation products, here the scalar h and pseudo-scalar a. The only decay final states relevant to antiproton production are actually quark pairs energetic enough to produce a pair of protonantiproton. This means that antiprotons can be created each time a light (pseudo-)scalar boson is heavier than twice the proton mass, roughly, which sets a production threshold to m φ 2 GeV.
The first step to get the antiproton spectrum is easy: before the decays of our light (pseudo-)scalars, we have a mere twoin-two-body annihilation so all the 4-momenta are set by kinematics. Assuming that dark matter particles annihilate at rest in the Galactic halo frame, the energy of resulting particle 1 is and the norm of its momentum is where , which are enough to change from the halo frame to the rest frame of particle 1: (γ, β)=(E 1 /m 1 , k 1 /E 1 ). In the rest frame of particle 1, the quarks it decays into have energy E * q = m 1 /2 and momentum |k * q | = m 2 1 /4 − m 2 q . One finally gets the energy of the quarks in the halo frame from a mere Lorentz transform: Thus, the energies of the quarks and anti-quarks coming from the decay of particle 1 are evenly distributed between E 1 2 −E and E 1 2 + E. One finally gets the probability of having an antiproton of energy Ep from a quark of energy E q f (E q , Ep) thanks to the PYTHIA 4 package [86]. Eventually, the antiproton spectrum obtained after decay of the annihilation products is: where the first sum is done over the annihilation products (1,2)=(a, h), and the second sum is over all those quark flavors for which 2 m q m i is satisfied. The branching ratios Br i,q can be derived from the couplings to SM fermions given in Eq. (5) and the decay widths given in Eq. (6). For the quark masses, we use the pole masses reported in e.g. [87]. It is clear that the decay into the down-type b quark-antiquark pairs will dominate whenever permitted. When m i < 2 m b , the dominant decay channel becomes φ i → τ + τ − because of the more favorable dependence into tan β with respect to the uptype quark c. Therefore, the antiproton production is reduced in the range 2 m τ m i 2 m b , but not completely suppressed, since d and s quark contributions remain significant. We have calculated the decay branching ratios explicitly in Fig. 5 assuming c ϕ = c θ = 0.9 and tan β = 10 (50) on the left (right) panel. Those associated with the third, second and first generation fermions are reported with solid, dashed, dotted lines, respectively. These curves are relevant to h and a, but we recall that we have only constrained the mixing of h so far. For simplicity, we take the same mixing angles for a, which can be considered as conservative for the antiproton production -for 4 For this work we made use of version 6.4.24 with CDF tune A. example, a slight increase of the up Higgs doublet content of a would reduce the production of τ leptons.
Finally, we also consider the possibility that h decays into 2 a, a channel that is open when m h ≥ 2 m a . This channel is relevant when the coupling λ aah /m h C φqq (see Eq. 6). Given the range we have considered for λ (see Eq. 23), this barely happens in our scan. Still, we note that if this channel is open and dominant, then the upper limit on the mixing parameter s ϕ can be relaxed (see Sect. 2.3); for simplicity, we still do not relax it.
An illustration of the spectrum determination is presented in Fig. 6, where we have used different combinations of the masses m χ , m a and m h and applied the Lorentz boost procedure described above -we did not include h → a a. We have considered five different cases, with m χ = 12 GeV for each. We first see that each time m a , m h > 2 GeV, the antiproton production is rather efficient and peaks with almost the same amplitude, within a factor of ∼ 3, with a sharp decrease. The brown curve (third from top at peak energy) shows the case m χ = m a = m h , which means that a and h both mostly decay at rest into b quarks. Nevertheless, these quarks are created with a relatively small momentum (E b = m a,h /2 = 6 GeV), which reduces the proton and antiproton production efficiency. The violet (upper) curve shows a more favorable situation in which m a = 12 GeV and m h = 4 GeV, i.e. both a and h can decay into quarks with momenta large enough to fragment into antiprotons efficiently -the center of mass energy is 2 m χ , and the production of τ is still not strongly dominant for the h decay because very close to threshold (see Eq. 6 and Fig. 5). The pink curve (second from top) characterizes an unexpected situation where, although h and a decays dominantly into τ leptons, the remaining ∼ 30% going into s and d quarks are very efficient at producing antiprotons thanks to their large momenta. The orange curve (fourth from top) shows a less favorable case for antiproton production in which a is not massive enough, and the b quarks coming from the h decay have low momenta, which translates into a factor of ∼ 2 suppression with respect to the previous case. Finally, the blue (lower) curve corresponds to an even worse case, too small m a and m h = 10 GeV, where b quarks are inefficient at producing antiprotons in addition to threshold effects (the decay into τ comes into playsee Fig. 5).
These spectral differences that depend on the mass combination demonstrate that there are non-trivial ways to suppress the antiproton production, e.g. producing b quarks at threshold with low momenta. Nevertheless, we will further show that each time a or h can decay into antiprotons efficiently, the antiproton flux is always in tension with the current data, provided the s-wave contribution to the annihilation cross section given in Eq. (12) is dominant at freeze out.

Dark matter halo
To predict the antiproton flux at the Earth, we need to specify the dark matter distribution in the Galaxy. Many conventional dark matter halo functions and related parameters have been  widely used in the literature so far, but we aim at employing a distribution as constrained as possible, both theoretically and observationally. From the theoretical point of view, there seems to be a consensus in using an Navarro-Frenk-White (NFW) profile [88,89] or an Einasto profile [90,91], which are found to provide good fits to the dark matter distributions in cosmological N-body simulations (see e.g. [92]). The former is mildly cuspy in the center, behaving like r −1 , while the latter asymptotically flattens towards the center. Nevertheless, we stress that the very central parts of the Galaxy could differ from what most of dark matter-only N-body simulations have predicted because baryons dominate the gravitational potential there, which might modify the central behavior of the distribution [93]. From the observational point of view, observations of low surface brightness galaxies also seem to indicate cores instead of cusps in these objects [94], showing that profiles are likely not universal in the centers. Still, in contrast to gamma-rays, the antiproton signal is well known to have almost no dependence in the distribution in the Galactic center, provided it is less cuspy than r −3/2 [69,79,32]: this is due to diffusion effects. Much more important is the local dark matter density, which roughly sets the antiproton flux. It is therefore of stronger relevance to use a model dynamically constrained on kinematic data.
Recently, complementary efforts were undertaken to constrain the local density of dark matter ρ [84,85,95], the results of which all point to ρ 0.4 GeV/cm 3 . Here, we will use the spherically symmetric NFW profile as constrained in [95], which we recall hereafter: the variable r denoting the galactocentric radius, and setting the parameters to r = 8.28 ± 0.16 kpc, ρ = 0.40 ± 0.04 GeV/cm 3 and r s = 18 ± 4.3 kpc. We emphasize that the antiproton flux will be sensitive mostly to variations in ρ , not in the other parameters.
In principle, consistency imposes to consider dark matter substructures [96,97], which can increase the antiproton flux at low energy [98,99]. While the amplification is expected to be small [99], we stress that even a factor of 2 could have important consequences in delineating the antiproton constraints, as will be shown later. Here, we will not consider substructures, which can therefore be considered as a conservative approach.

Antiproton transport and flux at the Earth
Once antiprotons are injected in the Galaxy, they diffuse on magnetic inhomogeneities and experience different additional processes. At low energy, say below a few GeV, they are sensitive to Galactic winds which translate into convection upwards and downwards the Galactic disk and into adiabatic energy losses; they are also sensitive to nuclear interactions with the interstellar gas, which may destroy them and/or decrease their energy; finally, they can be reaccelerated because of the proper motion of the magnetic scatterers. Higher energy antiprotons are roughly sensitive to spatial diffusion only. The general steady state equation that drives their evolution in phase-space reads: where N = N( x, p) = dn/d p = β dn/dE is their differential density, K xx and K pp are respectively the diffusion coefficients in space and momentum, V c is the convection velocity and Γ s = Γ s ( x, p) the spallation rate. This equation can be solved numerically [100] or semi-analytically [101,102,103]. In this study, we have used the semi-analytical method sketched in e.g. [103,69,79,98], to which we refer the reader for more details. For the transport parameters, we have utilized the bestfit model found in [103] and still in excellent agreement with more recent analyses (e.g. [83]). This transport model was dubbed med (for median) configuration in [69]. It is a slab transport model, with a cylindrical diffusion volume of radius R = 20 kpc, half-height L = 4 kpc, in which the diffusion coefficient obeys K xx (E) = β K 0 (R/1 GV) δ (R = p/q is the cosmicray rigidity). The main parameters are reported in Tab. 1. The source term is characterized by the specific dark matter model configuration: where σv = a ah is the annihilation cross section in the Galaxy given in Eq. (12), and dN dec /dE is the antiproton spectrum arising from the (pseudo-)scalar decay and determined in Sect. 4.1.
Assuming a Green function G( x, p ← x s , p s ) solution to Eq. (29), the antiproton flux at the Earth is given by the fol-lowing expression: In practice, we use the Bessel series method presented in [69,79], but we have checked with a couple of configurations that results obtained from Green function methods [79,98] are equivalent.

Results and discussion
We have performed a scan over the masses and couplings associated with Eq. (1) such that we keep only those configurations falling in the so-called CoGeNT region (see Sect. 3). For these configurations, we have computed the annihilation cross section summing over the three diagrams of Fig. 2 for all possible final states, i.e. a a, h h and a h as soon as kinematically allowed -we recall that we have imposed 2 m χ ≥ m a +m h . Then, we have derived the relic density assuming different QCD phase transition temperatures, and computed the antiproton fluxes for models leading to the correct relic abundance, i.e. 0.0941 ≤ Ω χ h 2 ≤ 0.1277 in this study (see Sect. 2.4). The flux is computed assuming the transport parameters of Tab. 1, and a solar modulation is subsequently accounted for in the force-field approximation [104,105], with a Fisk potential of Φ = 600 MV, consistent with the average solar activity during the PAMELA data taking [73]. For the secondary antiproton flux, we use the predictions obtained in [78] with the same transport parameters. It is represented in the left panel of Fig. 7 together with a prediction using another parameter set characterized by a larger diffusion halo with L = 15 kpc. Both predictions are consistent with the current PAMELA data [73], while the latter transport model would lead to a significantly larger primary antiproton flux because integrating over much more dark matter annihilation products. It turns out that large diffusion halo models, i.e. with L 5 kpc, are favored by several recent studies on Galactic cosmic-ray transport [83,81], while small halo models with L 2 are strongly disfavored from complementary observations [101,106,32] and are to be considered as extreme cases -we complete this discussion in Sect. 4.3.2. Therefore, we will stick to our median transport model which can be thought of as rather conservative.
We finally confront the predicted antiproton fluxes with the PAMELA data 5 . An error box is defined as having the energy bin width and the 1-σ error bar half-height. A model is considered in tension with the data when the predicted flux (adding up the primary and secondary contributions) is above an error box without crossing it 6 . A few examples are shown in the right panel of Fig. 7, where the same injected spectra as in Fig. 6 have been used an compared to the case of an annihilation into bb. A conventional annihilation cross section of σv = 3 × 10 −26 cm 3 /s was assumed, for which all models appear in excess with respect to the data. This clearly illustrates  Fig. 6; they are compared to a bb spectrum similar to the one used in [32]. the rejection power of the antiproton analysis. In the following we present the results we obtained from scans that cover the CoGeNT region, as detailed in Sect. 3. Fig. 4 shows the results of two scans we performed to cover the CoGeNT region in the plane m χ − σ SI . Dark red points are those which fulfil the relic density constraint, and the light orange points are those in excess with respect to the antiproton data. Both populations can roughly cover the whole region, for small (upper panels) as well as for large values (lower panels) of tan β, independently of T QCD . Still, taking small values of tan β results in more cosmologically relevant configurations in average, as it can be observed. Indeed, tan β = 10 implies a coupling c χh a factor of ∼ 10 larger than in the case of tan β = 50, for a given spin-independent cross section (see Eq. 21), which provides more room to get the correct relic density only from the annihilation into h h. Modifying the QCD phase transition temperature leaves these statements valid (see left/right panels). More generally, Fig. 4 demonstrates that we cannot draw a generic antiproton exclusion contour in this plane, since many configurations pass both the cosmological and the antiproton tests. Nevertheless, it shows that the antiproton constraints are relevant everywhere in this plane, which is by itself an important result. In the following, we provide details about the nature of those points which are (or not) in excess with respect to the antiproton data.

CoGeNT region scan analysis
In Fig. 8, we replot the scan results of Fig. 4 in the plane σv 0 -m χ . Note that σv 0 refers to the annihilation cross section relevant to indirect detection, i.e. the velocity-independent contribution given in Eq. (12), while the full velocity-dependent annihilation cross section was used in the relic density calculation. We remark that the relic density constraint obviously translates into an upper bound σv 0 5(8) × 10 −26 cm 3 /s for T QCD = 150(400) MeV, but that the cosmologically viable points spread down to arbitrarily small values of σv 0 that correspond to cases where the P-wave is dominant, i.e. where annihilation proceeds mostly into h h or a a at freeze out. This demonstrates that the kinematic condition 2 m χ ≥ m a + m h does not generically ensure that annihilation into a h dominates. Indeed, suppressing the coupling between the singlino and the pseudo-scalar a is enough to suppress the S-wave, as shown in Eq. (12). We will come back to this later. The population located close to the upper cosmological bound is featured by those configurations for which annihilation is saturated by the S-wave, i.e. annihilation into a h. This is precisely where the antiproton constraint is the strongest. There is still a significant population that remains unaffected by the antiproton constraints in this region. This actually corresponds to cases where the decays of a and h do not produce antiprotons efficiently (see discussion in Sect. 4.1). Eventually, the impact of changing T QCD is a bit harder to figure out from our scan results. Typically, considering large values of T QCD implies a larger cross section at freeze out over a broader WIMP mass range, i.e. when m χ /x f T QCD (see Fig. 1). This means that going from T QCD = 150 to 400 MeV, we increase the required annihilation cross section by a factor of ∼ 2 in the singlino mass range ∼4-12 GeV (see Fig. 1).
A way to identify the impact of the S-wave content of the annihilation cross section on the antiproton constraint strength consists in determining the S-wave to P-wave ratio, as shown in Fig. 9. Each panel illustrates the ratio a ah /b hh versus the ratio a ah /b aa for different combinations of tan β (10/50 on top/bottom panels) and T QCD (150/400 MeV on left/right panels). The additional factor of 20 applied to the ratios comes from that we where a ah , b aa , and b hh are defined in Eqs. (12,13,14), respectively. In each panel of Fig. 9, going from bottom to top means decreasing the relative contribution of the annihilation into h h, while going from left to right means decreasing the relative contribution of the annihilation into a a. The top left part of each panel, where annihilation into a a is always dominant with respect to other channels, is barely populated. In contrast, there is a dense wake of points going diagonally from the bottom right part to the bottom left part, which is charac-teristic of an annihilation driven by h h production; indeed, as long as SP hh < SP aa < 1, this channel overtops the others 7 . This implies that annihilation into a a is always subdominant in our sample. This is partly a consequence of requiring the Co-GeNT region to be fulfilled, which constrains the coupling c χh to be large, while c χa is left free; besides, we see that going to larger values of tan β slightly alleviates this effect, as expected (see Eq. (22)). Another cause is less trivial: increasing c χa does not only increase annihilation into a a, but also that into a h (see Eq. 12); the latter is always kinematically allowed here, but not necessarily the former. Fig. 9 further lets us distinguish a particular area, made Figure 9: S-wave to P-wave ratio, a ah /b aa versus a ah /b hh , as calculated from Eqs. (12,13,14).
shaded in each panel, that corresponds to S P aa > 1 and S P hh > 1. This area is the region where the S-wave saturates the full annihilation cross section at freeze out, and is therefore fully fixed by the relic density (optimal case for indirect detection). Naturally, the largest fraction of the models inducing an excess in the antiproton flux (light orange points) lies in this region. Nevertheless, there are still some configurations unaffected by the antiproton constraints there, which are actually those for which a and h are too light to decay into antiprotons (through quark fragmentation), or with masses 10 GeV for which the decay into b quarks occurs at threshold and is is then inefficient at producing antiprotons (∼ [9.5,12] GeV, see Sect. 4.1). In the latter case, m χ usually takes rather large values because of our kinematical condition, which also decreases the antiproton flux (scaling like 1/m 2 χ ). We note that the density of unconstrained configurations increases with tan β, which comes from three ef-fects. First, the up/down asymmetry in the branching ratio is strongly dependent in tan β, and taking a large tan β increases the relative decay into τ leptons compared to that into c and s quarks (see Eq. 6 and Fig. 5). Second, a larger value of tan β induces an increase in m h for a given spin-independent cross section (see Eq. 22), and either the density of points in the mass range m h ∼10-12 GeV increases accordingly and m χ is pushed towards larger values. Third, a larger value of tan β may also translate into a smaller χ χ h coupling, leaving more room to the h a a coupling to set the relic density from the s-channel (see Eq. 12 and the left diagram of Fig. 2). The consequence is that h → a a opens up when masses permit, which may further reduce the antiproton production efficiency.
In any case, the plots of Fig. 9 show that the antiproton constraints can be generically very strong whenever S P aa 1 and S P hh 1, as expected -one can safely use Eqs. (12,13,14) to account for this constraint.
In Fig. 10, which depicts the indirect detection annihilation cross section σv 0 as a function of the singlino mass m χ , we report the same samples as in Fig. 8 after applying the cut S P tot > 1 (see Eq. 32). We clearly see that this cut selects all points for which the relic density completely sets σv 0 , a large fraction of which turns out to lead to an antiproton excess. The remaining unconstrained cases are mostly characterized either by m a 2 GeV, for which a cannot decay into antiprotons, or the condition m χ 10 GeV, mostly due to the antiproton flux dependence on 1/m 2 χ (see Eq. 30). This appears more clearly in the plane m h -m a shown in Fig. 11, where our results are presented with and without the cut S P tot > 1 -a rectangle is drawn in each panel to emphasize the most constrained mass range. We see that once the cut is applied, all models exhibit an excess with respect to the antiproton data but those with m a 2 GeV, and (m a + m h ) 20 GeV -the latter condition is equivalent to requiring m χ 10 GeV.
The dependence on tan β and T QCD arises in all Figs. 8, 9, 10 and 11. Small tan β values are generically more constrained by the antiproton data because they lead to relatively small values of m h for a given spin-independent cross section (see Eq. 21), thus favoring smaller values of m χ which increases the antiproton flux. In this regime, the a/h decay into τ leptons is also slightly less significant than in the large tan β regime. As to the impact of T QCD , things are much simpler. The cross section needs to be relatively larger for m χ 20 T QCD , which makes the indirect detection constraints more stringent for large values of T QCD .
Coming back to the full results shown in Fig. 8 in the light of the previous discussion, we may summarize our results with a generic upper limit on the annihilation cross section relevant to for This is not an accurate limit, but rather a way to avoid regions of the parameter space that could lead to an antiproton excess. This is illustrated in Fig. 12, where we used the worst situation for the antiproton constraints, tan β = 50 and T QCD = 150 MeV. We see that a looser cut of S P tot > 0.3 still leads to an excess in the specified mass range. Finally, we show in Fig. 13 how the antiproton constraint may translate into a constraint on the spin-independent cross section (see Eq. 21), comparing full samples (top row panels) to the same ones but after applying the cut S P tot > 1 (middle row panels). In each top panel, the region where the relic density is set by χ χ → h h, P-wave case irrelevant to indirect detection, covers almost the full area but the extreme right-hand side, where χ χ → a h becomes dominant. The (anti)correlation between σv 0 ∝ [c χχa /m χ ] 2 and σ SI ∝ [c χh /(m 2 h c β )] 2 appears precisely there, which is more explicit in the middle-row panels, where the cut S P tot > 1 was applied. From these plots, one can extract the following approximate upper limit on the direct detection phase space: In the bottom row panels of Fig. 13, we report a more explicit translation of our antiproton constraints in terms of spinindependent cross section, which can be directly compared with Fig. 4. We have applied the cut S P tot > 1. We see that antiprotons provide quite strong constraints for m χ 10 GeV in this regime, independently of tan β and T QCD . This nicely illustrates why the antiproton signal must be considered carefully in the singlino-like phenomenology.

Uncertainties and potential detection prospects for
PAMELA and AMS-02 There are different sources of theoretical uncertainties that might affect our predictions. Among important effects which are related to completely different domains, some uncertainties may be connected to the QCD phase transition in the early uni- verse (relic density), the strangeness content of nucleons (spinindependent cross section), the solar modulation modeling and the cosmic-ray transport parameters (antiproton flux).
Uncertainties in the QCD phase transition have been shown to lead to a factor of 2 uncertainty in the annihilation cross section at freeze out in the WIMP mass range ∼ 4-12 GeV, with larger transition temperatures inducing an annihilation cross section for m χ /x f T QCD larger by a factor of 2 than the one obtained for m χ /x f T QCD , due to current constraints on the dark matter abundance (see Fig. 1). We tried to account for this uncertainty by using both a sharp QCD phase transition and two different temperatures.
Regarding uncertainties in the hadronic matrix elements, we refer the reader to Refs. [60,61,62,63,64] for more details. The consequences on our analysis are twofold, since these uncertainties can lead to select either higher (larger strangeness fraction) or smaller (smaller strangeness fraction) values of C χh /(c β m 2 h ) -see Eq. (22). In the former case, then annihilation into h h through the t-channel would be increased, which would likely make the population obeying SP hh < SP aa < 1 more prominent, thereby leading to a P-wave dominated annihilation cross section. This would further squeeze the configuration range available to indirect detection. Nevertheless, we point out that requiring a correct relic density poses significant constraints then, since the annihilation rate may increase too much and lead to a too small dark matter abundance. On that account, it may become difficult to find cosmologically viable configurations, especially at small tan βthis is illustrated in the top panels of Fig. 13, where the upper regions are shown less favored by the relic density constraint. In contrast, the latter case, namely selecting smaller values of C χh /(c β m 2 h ), is more favorable to indirect detection, since a correct relic density can be achieved by means of a larger coupling C χa , increasing the annihilation into a h; so that the antiproton constraints get stronger.
Coming back to astrophysics, the low energy cosmic-ray flux below a few GeV is sensitive to the solar activity. In our calculation, the solar modulation is taken into account by using the force-field approximation [104,105], which is known to provide a rather satisfying understanding and fit, though missing some details, of the impact of solar activity on cosmic rays [107,108,109]. We refer the reader to Ref. [109] for a more complete review about the current state-of-the-art. Most of the solar modulation models are based on numerical or analytical solution to the Parker transport equation [110], the forcefield approximation providing an analytical solution by assuming perfect spherical symmetry and a static electromagnetic field. More complete models lead to charge sign effects depending on the solar magnetic polarity, which reverts at each 11-year solar cycle, but they are expected to be small, at the ∼10-20% level in terms of flux at the Earth. This roughly sizes the theoretical uncertainties affecting the force-field approximation. This has to be compared with the error bar associated with the lowest energy PAMELA data points, typically 1.5 GeV, which turn out to be the most constraining in our analysis (see Fig. 7)the observational error is ∼ 40% in this energy range, which is twice more important than the solar modulation error. We have tuned our Fisk potential to a low/intermediate value of φ = 600 MV, a rather conservative assumption with respect to the solar activity monitored during the PAMELA data taking [111]. This means that the error boxes that we have considered to derive our constraints safely encompass the theoretical uncertainty affecting the solar modulation modeling. In Fig. 14, we show the ratio of the total flux to the secondary flux, which aims at illustrating the spectral break due to additional primaries for differ-ent amplitudes characterized by the annihilation cross section σv = 10 −26 cm 3 /s (10 −27 cm 3 /s) on top (bottom) panels -and for different solar activities -φ = 600 MV (900 MV) on left (right) panels. We see that going from weak to strong solar activity can decrease a predicted excess down to ∼ 40%. The PAMELA data, which were taken close to a solar minimum, are therefore very well suited for to derive low energy limits.
As regards the uncertainties affecting the cosmic-ray transport parameters, we stress that one of the most sensitive parameter for dark matter-induced cosmic rays is the diffusion halo size L, which characterizes the volume which the annihilation products are integrated over and which turns out to be strongly correlated with the diffusion coefficient normalization K 0 -K 0 /L ≈ cst from stable secondary-to-primary nuclei ra-tio constraints like boron-to-carbon (e.g. [112]). Increasing L leads to integrate more annihilation products and thus increase the dark matter-induced primary flux [69,98,113]. In this paper, we have used L = 4 kpc (see Tab. 1). We emphasize that the most recent works on fitting the transport parameters tend to favor a larger diffusion zone, with L ∼ 5-10 kpc, which makes our propagation model rather conservative. These studies consider the additional input of the secondary radioactive cosmic ray species which allow to break the degeneracy between L and K 0 [83,81]. Another unexpected problem in small-L diffusion zone models is that they systematically lead to a secondary positron flux in significant excess with respect to the data at low energy [114,115]. Few other arguments in favor of larger values of L are recalled in [32] -see complementary arguments in a recent study on the diffuse Galactic radio emission [116]. Our values for K 0 and L can therefore be considered as quite motivated and rather conservative.
Finally, it is interesting to make some prospects for the end of the PAMELA campaign or the new operating space experiment AMS-02 8 , which was successfully launched in May 2011 and which is supposed to take data for 10 years. For the latter, it is not yet clear whether the low energy systematic errors can be reduced with respect to PAMELA, but statistics will for sure be much better. Fig. 14 shows how a finer PAMELA/AMS-02 data analysis could improve our current limits or lead to the detection of those singlino configurations which still induce cosmic antiprotons, but with a flux too weak to be detected or constrained. This typically corresponds to an annihilation cross section σv 5 × 10 −27 cm 3 /s. We see in the bottom panels that a precision of ∼ 5% makes the data quite relevant to such a low annihilation cross section, even in the case of strong solar activity. Therefore, if the PAMELA or AMS-02 collaborations manage to reach such a low precision, we expect much tighter constraints on light dark matter models annihilating into antiprotons in the future, even if the solar activity is also expected to increase in the coming years. In particular, much deeper phase-space regions with S P tot 0.1 (see Eq. 32) could be probed.

4.3.3.
Connecting this effective approach to more constrained particle phenomenologies The question may arise of how to connect this effective approach to existing and more constrained particle phenomenologies beyond the SM.
First, an interesting and very general link can be made by focusing on the kinematics of the final states. Indeed, we have considered that annihilation proceeds into two bodies that may further decay into quarks, leading to the production of antiprotons. We have shown that each time the annihilation cross section is S-wave dominated, then this kinematical configuration may lead to an antiproton excess quite generically. Our effective approach is therefore fully relevant to cases for which such two-exotic-body production is predicted, as in e.g. Ref. [117], where the so-called DAMA/CoGeNT signal is explained from a model incorporating two additional light gauge bosons that dark matter can annihilate into. Since one of these gauge fields has a mass larger than 2 GeV, it can in principle decay into an antiproton. This model can therefore suffer the antiproton constraints. Our constraints may also be important for completely different models, like in the recent study performed in [118] of a light right-handed sneutrino dark matter that can arise in the frame of the NMSSM; annihilation can then indeed go into light pseudo-scalars without being velocity suppressed, while this case is not necessarily adapted to fit the CoGeNT region.
Given our effective Lagrangian of Eq. (1), a natural link can be made with singlet extensions of the MSSM, which was the original motivation for this work. This includes general extensions as well as more constrained ones, like for instance the NMSSM, in which the dark matter candidate is a singlinodominated neutralino. An example of a model of singlino dark matter can be found in Ref. [23], where the specific annihilation into light (pseudo)scalar Higgs bosons is considered. Their case (3), which corresponds to 2m χ > m a + m h , can be directly compared with our analysis. In this specific configuration, they find |C χh | ≈ |C χa | |λ i jk |, with |λ hhh | |λ aah |. Having similar couplings between χχa and χχh implies that they both have to be large to explain the DAMA/CoGeNT signal, which also implies that the S-wave contribution of the t-channel annihilation diagram into a h is always significant. The fact that |λ hhh | |λ aah | induces a competition with the s-channel annihilation into h h (P-wave), but more suppressed than the t-channel annihilation into a h (velocity suppression, and a squared propagator factor of ∼ 1/(4m 2 χ −m 2 h ) 2 as compared with ∼ 1/m 4 χ for the t-channel). The antiproton limit must therefore be regarded in this case.
For comparisons with the NMSSM, we may use Refs. [22,119,120,121]. Although the reason is not clear from the mass matrix elements only [25], it appears to be difficult to generically find viable NMSSM configurations obeying 2m χ > m a + m h . Bosons a and h can be found very light, but not simultaneously [119,120,121] -an exception can be found in Ref. [22], but the authors get an annihilation going dominantly into ff . The above cited studies still find complementary viable regions in the NMSSM parameter space. When the relic density is set from the s-channel resonant exchange of a light h, then the annihilation is P-wave dominated and indirect detection fails to yield interesting bounds. If, in contrast, the main annihilation final state consists in an S-wave production of quarks (even if not with a branching ratio of 1), then an antiproton excess very likely follows in most of casesthis may also happen for a right-handed sneutrino dark matter within the NMSSM [118]. This situation was already discussed in Ref. [32], where the antiproton limit was shown to be very sharp.

Conclusion
In this paper, we have designed a generic effective Lagrangian to survey the light singlino dark matter phenomenology, as soon as annihilation mostly proceeds into light singletdominated scalar h and/or pseudo-scalar a Higgs bosons (see Sect. 2). Within this general framework, we have computed the annihilation cross section (Sect. 2.5), the relic density (Sect. 2.4), the spin-independent nucleon-singlino cross section (Sect. 3), the antiproton spectra arising from the decays of these light (pseudo-)scalar bosons (Sect. 4.1), and the subsequent antiproton flux at the Earth (Sect. 4.2.2). We have focused on a specific mass configuration that guarantees the existence of an S-wave contribution to the annihilation cross section, ensuring the relevance of indirect detection, i.e. 2 m χ ≥ m a + m h . We then have selected the mass and coupling parameter space to get a spin-independent cross section range encompassing the so-called CoGeNT region, and drawn several samples of ∼ 10 5 models for two different values of tan β, 10 and 50, for which we have calculated the relic density. Then, we have computed the antiproton fluxes associated with these configurations and compared them to the PAMELA data (Sect. 4.3).
We emphasize that there are still theoretical uncertainties affecting the estimate of the annihilation cross section at freeze out for ∼ 10 GeV mass WIMPs, up to a factor of ∼ 2 depending on whether WIMPs decouple after or before the QCD phase transition. This is due to the rapid change in the relativistic degrees of freedom that occurs at that time and acts as a weight on the annihilation cross section. We have taken this uncertainty into account by considering extreme cases, as illustrated in Fig. 1.
Our results demonstrate (Sect. 4.3) that the antiproton constraint is generically relevant to the whole CoGeNT region (see Fig. 4), though there are some simple ways to escape it. The most trivial way to suppress the S-wave contribution to the annihilation cross section is to suppress the coupling between the singlino and the light pseudo-scalar c χa , as shown by Eq. (12). Another way is to have the light scalar and pseudoscalar Higgs bosons so light that they cannot decay into antiprotons, or cannot generate antiprotons energetic enough to be observed. Anyway, we have shown that as soon as the S-wave is significant at freeze out, the antiproton signal is in most of the cases in excess with respect to the PAMELA data (see Fig. 8, which translates into Fig. 10 after imposing the S-to-P wave ratio to be larger than 1). Depending on the singlino/scalars mass configuration, we have derived a quite generic upper limit on the annihilation cross section lying in the range σv 0 10 −26 cm 3 /s, valid for the singlino mass range ∼6-15 GeV. This limit is very competitive with respect to those obtained from gamma-rays [30], and complementary to limits coming from high-energy neutrinos from the Sun [31]. We recall that the antiproton bounds are even more stringent for direct annihilation into quarks (see Fig. 7 and [32]).
Uncertainties affecting our analysis have been discussed in Sect. 4.3.2, where we have shown that our procedure is to be considered as rather conservative. In addition to the points raised there, we remind that our antiproton flux calculation is poorly sensitive to the details of the dark matter distribution far away from the Earth, in contrast with the calculation of the gamma-ray flux. The main parameter which sets the flux amplitude is the local density ρ , so one can easily rescale our results and limits by a mere factor of (ρ new /ρ ) 2 . Still, we recall that we have considered a dynamically constrained profile, and that we did not include any substructure. Although their survival against tidal effects and encounters with baryonic systems (stars and/or Galactic disk) still needs to be clarified, the survival fraction is expected to be larger and larger for smaller and smaller substructures, which are more concentrated [122]. Self-consistency would then require to include their effect in our calculation [97], which turns out to be a slight increase of the antiproton flux at low energy [98,99]. This would therefore make our present constraints even more stringent.
Finally, we have shown that provided systematic errors can be reduced down to ∼ 5%, PAMELA, from an improved analysis, or AMS-02 could have the capability of detecting a still unseen signal characterized by a spectral break (see Fig. 14) which would hide just within the current PAMELA error boxes we have considered in this paper. A null detection with better errors would in any case significantly ameliorate the limits derived in this paper.
We also introduce the functions g(s) and h(s) such that in the center of mass system: where f i j (s) is further defined in Eq. (A.5), and θ is the angle between p 1 and k i in the center of mass frame. The differential annihilation cross section is given by: , where g χ is the number of degrees of freedom of particle χ (g χ = 2 here), S i j ≡ 1 + δ i j is a symmetry factor, and where the integrated Lorentz invariant phase space factor reads: Note that z ≡ cos(θ).
We are armed enough to derive the cross section in the most general case, i.e. we can take particles φ i, j,k having all a scalar 9 We stress that this general writing introduces non-physical degrees of freedom when going to the total squared amplitude, notably for the interference terms (mixing of two different final states, e.g. scalar-scalar with pseudoscalar-scalar). One has to remember that interaction terms involving c χi ×c χi are forbidden -this is clear from the Lagrangian of Eq. (1) -and can therefore be dropped. As a more general rule, one can replace c χi ×c χ j −→ c χi ×c χ j (1−δ i j ). Despite this caution, which is relevant when trying to give physical interpretations to the general-case result, such a general approach proves powerful since a single expression can be used for different annihilation final states; one has just to switch the appropriate couplings on/off to get the selected one, which automatically removes the unphysical terms. and a pseudo-scalar component; we will take care of removing the non-physical degrees of freedom at the end 9 . Denoting σ β = σ 12 β 12 and using g = g(s) and h = h(s), we find: with α mt → (−1) m α mt ; and index i ↔ index j for odd m Parameter K, which carries dimensions of inverse squared mass, is defined as The coefficients α are associated with the numerators of terms involving the t-or u-channel, which all exhibit an angular dependence. This angular dependence can be expressed as a polynomial function of z reading P x (z) = n α nx [h(s) × z] n . These coefficients are explicitly given below: where we have used: 12c) The null-velocity limit of the annihilation cross section corresponds to requiring s → 4 m χ 2 . This implies that function h(s) → 0, which a posteriori justifies the form of Eq. (A.7), from which the full limit can be readily calculated. In terms of the above α coefficients (we now denote themᾱ to make explicit that they also have to be taken in the limit s → 4 m χ 2 ), we obtain the following suitable expressions: The above expressions are the ones which are relevant to indirect detection. We note that for the s-channel, only the pseudoscalar exchange gives a non-zero contribution, as expected. It is also easy to check from these expressions that only annihilation into a scalar and a pseudo-scalar is non-zero in the null velocity limit.
To get the Taylor expansion terms of the thermally averaged cross section of Eq. (11), we refer the reader to Ref. [44]. Defining the Lorentz invariant function The zeroth order term a (often called S-wave, though it is only a part of it) is therefore fully determined by Eq. (A.13), while the first order term b requires a bit more of algebra. The latter case is left to the patience of the reader, though explicit results are given for the cases of χ χ → h h and χ χ → a a in Eqs. (13) and (14), respectively. For the indirect detection signal, we have used the analytical expression obtained for the Swave, while for the relic density calculation we have computed the thermal average numerically following Refs. [46,47]. We recall the full form of the thermally averaged annihilation cross section below: σv = d 3 p 1 2 E 1 (2 π) 3 d 3 p 2 2 E 2 (2 π) 3 f χ ( p 1 ) f χ ( p 2 ) W(s) , (A. 16) where function f is the WIMP phase-space distribution, usually taken as a Maxwell-Boltzmann distribution. The above expression can be further developed and simplified, and is finally found to be a temperature-dependent integral over simply two variables: s, the squared center-of-mass energy, and µ, the cosine of the angle between the annihilating particlessee [46,47] for more details. As shown above, the integral over µ is performed analytically in our case. In practice, we compute σv numerically from Eq. (A. 16), where function W(s) is tabulated up to s − 4 m χ 2 = 10 m χ , including all potential resonance and thresold effects. The Taylor expansion provided in Eq. (11) is only meant to clarify the discussion. It is indeed well known that such an expansion is a bad approximation in some cases [123]. Still, the S-wave expression given in Eq. (12) turns out to be quite accurate to compute the indirect detection signals, which is not surprising given that the WIMP velocity distribution in the Galaxy peaks around x = m χ /T ≈ 10 7 .

Appendix B. Relic density
We refer the reader to Refs. [46,47] for more details on the relic density calculation, which we have followed in this study. We still provide a few pieces of information below, mostly related to the effective degrees of freedom, which is an independent part, except for the last paragraph which gives some information about the numerical algorithm.
Assuming entropy conservation and a Hubble rate set by the energy density ρ as H 2 = 8πρ/3, the Boltzmann equation that drives the evolution of the dark matter fluid of density n = s Y in the early universe reads: where M p is the Planck mass, s = (ρ + p)/3 is the entropy density (p being the pressure). The dynamical variable is x = m χ /T , increasing when the universe expands. We further define the effective relativistic degrees of freedom g eff and h eff relevant to the energy density and the entropy density, respectively, such that: ρ = π 2 30 g eff (T ) T 4 (B.2) s = 2 π 2 45 h eff (T ) T 3 It is useful to recall the expressions of the energy density and the entropy density, for a fermion or a boson species labelled i: exp E T ∓ 1 These integrals can actually be evaluated quite easily from series involving modified Bessel functions of order n, K n , which quickly converge to the exact results such that each species i carries [124]: where x i = m i /T . Note that in the ultra-relativistic limit x i → 0, these expressions tend to g i and (7/8) g i , respectively for bosons and fermions, as expected. We have used the following property of Bessel functions [125]: +∞ y x 2 − y 2 α−1 e −µ x dx = 1 √ π 2 y µ 2 α−1 2 Γ(α) K 2 α−1 2 (µ y) .
Given these expressions, one can calculate the complete set of effective degrees of freedom by summing over all physical degrees of freedom at a given temperature. Before the QCD phase transition in the very early universe down to T ∼ M W , this amounts to all degrees of freedom of the standard model, including quarks and gluons, except the W ± , Z and Higgs bosons (the electroweak phase transition occurs earlier when T M W ). After the QCD phase transition, quarks and gluons are no longer free and remain confined into hadrons; then resonances have also to be included in addition to nucleons, an can provide extra-contribution depending on whether they are relativistic or not [126].
Note that in principle, one should also add those additional degrees of freedom coming from non-standard particles related to the particular dark matter model under scrutiny. In this paper, the light Higgs bosons and the singlino should therefore be added. Nevertheless, we can safely neglect them since they are also non-relativistic when singlinos decouple, around m χ /T ∼ 20 (we used m a 1 GeV and m h 4 GeV, and focused on 5 GeV m χ 15 GeV). Our results are illustrated in Fig. 1, where a sharp QCD phase transition has been considered. Such a transition implies that g goes to infinity because of the derivative of h eff in Eq. (B.2). We have suppressed this divergence based on that the actual transition is expected to be smoother [54]. We have checked that considering this effect or not has no significant impact on our results, at the percent level.
As regards the numerical method, we have implemented a Runge-Kutta-6 scheme to track the evolution of the comoving density Y featuring Eq. (B.1). To avoid starting the calculation too early (too small x), we first estimate the decoupling temperature by implicitly solving Eq. (B.1), demanding Y(x ) = (1+ )Y eq (x ) (with < 1/2, typically) -before decoupling, Y tracks its equilibrium value. We then start the calculation at x in x , up to after the decoupling. Indeed, to achieve an accuracy at the percent level, we need to solve Eq. (B.1) up to x end such that Y(x end ) > 10 Y eq (x end ), because of the term (Y 2 − Y 2 eq ). In practice, we use Y(x end ) = 30 Y eq (x end ). Then, we neglect Y eq and can therefore use the solution of Eq. (9), taking x > f = x end and performing a Simpson integral over x in logarithmic steps.