Fermion Singlet Dark Matter in a Pseudoscalar Dark Matter Portal

We explore a simple extension to the Standard Model containing two gauge singlets: a Dirac fermion and a real pseudoscalar. In some regions of the parameter space both singlets are stable without the necessity of additional symmetries, then becoming a possible two-component dark matter model. We study the relic abundance production via freeze-out, with the latter determined by annihilations, conversions and semi-annihilations. Experimental constraints from invisible Higgs decay, dark matter relic abundance and direct/indirect detection are studied. We found three viable regions of the parameter space, and the model is sensitive to indirect searches.


Introduction
Astrophysical evidence of dark matter (DM) has been accumulating for more than forty years now, but its fundamental nature remains unknown. From the particle physics points of view, different approaches have been carried out over the years to account for the elusive DM (for a review see [1]), and in particular, simple extensions to the SM containing gauge singlets look appealing for their simplicity, DM predictions and testable phenomenology [2][3][4][5][6][7][8]. Nowadays, those WIMP minimal extensions have been very constrained, especially by direct detection, motivating other possible alternatives. For instance, the interplay of a (pseudo)scalar and a fermion, both gauge singlets, open up the possibilities in many aspects: multi-component DM, new interaction channels, novel experimental signatures, small-scale structures, among others .
Keeping minimality, in this work we study a gauge singlet sector composed of a real pseudoscalar and a Dirac fermion. Depending on the coupling values and mass hierarchy between the singlets, the model admits a variety of DM productions with either one or the two singlets being stable. Two new interactions are present in the model: a Higgs portal and a dark sector coupling. The Higgs interaction is key because it regulates to what extent the dark sector is coupled to the SM, whereas the internal dark sector coupling only regulates the coupling between the two singlets. In our knowledge, we study for the first time the WIMP regime of this framework in which both couplings take sizable values such that both singlets were in thermal equilibrium with the SM bath in the early universe. Interestingly, in certain mass hierarchy between the two singlets, the stability of both fields is guaranteed by a parity symmetry without the necessity of introducing new ad hoc discrete symmetries. Models with this last feature or accidental symmetries have been studied in different DM context, such as Minimal Dark Matter [30], spontaneous symmetry breaking [31], two DM components [32], vector DM [33,34] and rank-two fields [35].
In the model under consideration, the DM relic abundance is triggered by annihilations, DM conversions [36] and semi-annihilations [37,38], showing remarkable features in some regions of the parameter space. Further, we constraint the model considering the measured relic abundance in the universe, Higgs invisible decay and direct/indirect detection. For the latter, we explore box-shaped gamma ray spectra [39,40], and we confront the available parameter space with Fermi-LAT data, CTA projections and AMS-02 bounds.
The paper is organized as follows. In section 2 we present the model and its theoretical constraints. In section 3 we explore the possible DM relic abundance mechanisms presents in the model, with a precise analysis of the two-component freeze-out scenario. In section 4 we review experimental constraints and the available parameter space along with indirect detection signals. In the last section we discuss and state our conclusions.

The model
The model adds to the SM two gauge singlets: one Dirac fermion ψ and a pseudo-scalar s. Under a parity transformation, the fields transform as ψ → γ 0 ψ and s → −s, giving rise to the following Lagrangian: where the scalar potential is given by with H being the Higgs doublet 1 . We make g ψ real under a fermion field redefinition such that the theory is CP-conserving. We assume that the singlet scalar does not acquire vacuum expectation value (vev), and after EWSB in the unitary gauge, GeV the Higgs vev, the scalar potential may be rewritten as with the mass of the scalars given by Here we consider m h = 125 GeV. The stability of the fermion is easily recognized due to the fact that it appears in pairs in the Lagrangian. In the scalar sector s appears in pairs, forbidding its decay. The linear term in s in 1, imply that as m s ≥ 2m ψ the pseudoscalar may decay into a pair of ψ, whereas as m s < 2m ψ the scalar singlet becomes stable at tree level and at all orders in perturbation theory. In the following we argument the latter fact. The possible decay of s at an arbitrary number of loops is represented in Fig. 1, with the decay of s followed by a singlet fermion closed loop and the black circle representing possible interactions between an arbitrary number of s and h. In the figure n , n and N simply represent the number of scalar lines. For simplicity, let us start assuming n = 0. If n is even, the resulting fermion trace will at most contain terms of the form µνρ... p µ p ν p ρ . . . , which vanishes exactly. If n is odd, the trace is different than zero, but there is no way to connect an odd number of outgoing s with an arbitrary number of N Higgs boson in the black bloop, due to the presence of the CP symmetry of the scalar potential. Now, if n = 0, the previous arguments still remains, because internal lines of s in the fermion loop only add an even number of both γ 5 and fermion propagators to the trace, only adding vanishing contributions without changing the final result. In consequence, from a perturbative point of view, the stability is guaranteed for the pseudoscalar singlet 2 . Finally, theoretical constraints put bounds on the free parameters of the model. The stability of the electroweak vacuum for s and h imposes that [45] with λ H = m 2 h /(2v 2 H ) 0.1. From perturbativity we set that |g ψ |, |λ s | < 4π to ensure that loop corrections are smaller than tree-level processes [45]. Unitarity constraints are less stringent than the upper limits based on our perturbativity criteria [5,46]. The signs of g ψ and λ hs are not relevant for the analysis in this work due to the fact that the relevant processes depend quadratically on them, therefore the theoretical constraints set 0 < λ hs < 4π and 0 < g ψ < 4π.

Parameter Space
Depending on the intensity of the couplings and the mass hierarchy between m ψ and m s , the model may present different DM scenarios with one or two stable particles. When λ hs takes very small values, i.e. λ hs ∼ 10 −12 − 10 −6 , the singlet sector never thermalize with the SM, then the DM production may occurs via freeze-in and/or dark freeze-out [42,47]. On the contrary, as the Higgs portal coupling takes sizable values, the singlet scalar enters into thermal equilibrium with the SM (this work). Based on the latter fact, two possible coupling regimes concerning g ψ may be present: • g ψ 10 −6 : The two singlets in the dark sector will interact feebly. If m s > 2m ψ , s becomes unstable, and ψ may be produced via freeze-in through the decay of s and 2 → 2 scattering processes. The green region in Fig. 2 shows this parameter space in terms of the mass hierarchy. For m s ≤ 2m ψ , s becomes stable and define a DM candidate identical to that of the Singlet Higgs Portal model (SHP) [5,8], since ψ does not interfere in the dynamic of the former due to their feebly interactions. The SHP model has been exhaustively studied previously, and displays a highly constrained parameter space around the EW scale.
• g ψ 10 −6 : s will bring ψ fast into the thermal equilibrium. From the relic abundance point of view, the only relevant case here is when both singlets are stable, i.e. m s < 2m ψ (orange and blue region in Fig. 2 ), otherwise ψ would not have any channel to annihilate, giving rise to an overabundance (see Feynman diagrams in Fig. 3). Based on the latter point, different type of interactions appear that determine the relic abundance of each singlet.
In this work, we focus on this last scenario, in which both singlets are stable. It is worth to mention that even when the first case with one DM candidate via freeze-in is a perfectly viable DM model, it presents a challenge phenomenology [48] (see however [49]). Recently, novel scenarios with inverse semi-production have been proposed with interesting dark matter production and phenomenology [50,51].

Boltzmann equation
As we pointed out before, the two-component scenario of our interest occurs for m s < 2m ψ and when the couplings (g ψ , λ hs ) are sizable. In this case, at temperatures higher than the individual masses of the singlet, both DM components are in thermal equilibrium with the SM. The departure of the equilibrium occurs once the temperature goes below the masses of the singlets, and three types of scattering participate in this process: annihilations, semi-annihilations, and dark matter conversions (Fig. 3). Based on [52], the evolution of the individual singlet abundances Y i ≡ n i /s, with i = ψ, s, as a function of the temperature x ≡ µ/T , with µ = m ψ m s /(m ψ + m s ), is given by where we have defined with X referring to a SM particle, σv the thermally averaged cross section, and the entropy density s and Hubble rate H in a radiation dominated universe given by where G is the Newton gravitational constant, and g * and g s * are the effective degrees of freedom contributing respectively to the energy and the entropy density at temperature T [53]. The equilibrium densities, Y i,e ≡ n i,e /s, are calculated using the Maxwell-Boltzmann distribution, whose number density is given by with g i the internal spin degrees of freedom, and K 2 is the modified Bessel function of the second kind. The DM density parameter is given by where Y i,0 is the yield of each DM component today, i.e. after freeze-out.
In the following we analyze the behavior of the relic abundance obtained by solving Eq. (7) making use of micrOMEGAs 5.2.7a [52], which calculates automatically all σv in Eq. 7. We pay special attention to those regions in which Ω s h 2 acquires smaller values than in the SHP [5] for a fixed (m s , λ hs ). This will be a necessary requirement in order to evade stringent direct detection bounds on the pseudoscalar DM from a few GeV to the TeV scale. Due to the subtle behavior of semiannihilations in the t-channel, we divide the analysis in two parts, depending on the relative hierarchy between m s and m h . Finally, for convenience, we introduce a new parameter ∆ ≡ 2m ψ − m s − m h , such that as ∆ ≥ 0 s-channel semi-annihilations are present in the relic abundance calculation. This is represented in blue in Fig. 2 .

m s < m h
In this region all the processes in Fig. 3 may be present but the semi-annihilation in the t-channel. Now, when m ψ is small enough, i.e. ∆ < 0, then the s-channel semi-annihilation is not present either, and only pseudoscalar annihilations and DM conversions are present. Normally in this regime, an overabundance occurs as m s > m ψ , since ψ does not have effective annihilation channels. As m ψ > m s , DM conversions of the type ψψ → ss become effective, decreasing the overabundance of Ω ψ . If m ψ is big enough such that ∆ ≥ 0, semi-annihilations in the s-channel become efficient, decreasing Ω ψ even more. These characteristics can be seen in Fig. 4 (above), where the abundances are depicted as a function of m ψ , for m s = 40 GeV, g ψ = (0.1, 1, 10) and λ hs = 1, with the green and red regions representing s unstable and DM overabundance, respectively. The effects of conversions at m ψ = m s becomes sharper as g ψ increases, since σ ψψss v ∼ g 4 ψ (see Appx. A). Equivalently, schannel semi-annihilations near ∆ = 0 pushes down Ω ψ , and since both conversions and s-channel semi-annihilations depends inversely on m 2 ψ , their effectiveness decreases with the growing of m ψ , then making Ω ψ to increase for higher values of m ψ . The changes of Ω s are stronger only for m s > m ψ , where ss → ψψ conversions are effective, and due to σ ssψψ v ∝ g 4 ψ , Ω s decreases notoriously for such high couplings g ψ shown in the Fig. 4 (above) right plot. We do not show the dependence of the relic abundances on λ hs , but the changes are minor due to the fact that this parameter affects only the semi-annihilation intensity and modulates Ω s .
In Fig. 4 (below ) we show the abundances as a function of g ψ , keeping λ hs = 1. The relic abundance behavior depends strongly on the mass hierarchy and the magnitude of g ψ . As it was previously discussed, for m ψ < m s (left plot) it is Ω ψ that dominates the relic, independently of g ψ . In this case, λ hs would only change the relative abundance of the scalar singlet, but keeping very high values for Ω ψ . Oppositely, in the cases m ψ > m s (middle and right plot), the relic hierarchy do depends on g ψ , showing the effectiveness of conversions and s-channel semi-annihilation, respectively, with a notorious fall of Ω ψ in each case.

m s > m h
In this case, the t-channel semi-annihilation ψ + s → ψ + h is present, and it may participates strongly in the determination of Y s for ∆ < 0. The effectiveness of this semi-annihilation on Y s becomes highly sharp, due to the fact that once both DM components decouple from the SM plasma, Y s follows in a good approximation assuming Y ψ ≈ constant and Y s,e ≈ 0. The solution of Eq. (12) gives an exponential suppression for Y s , highly sensitive to the mass difference between the singlets. As an example of this behavior, in Fig. 5 we show the evolution of the densities Y ψ (blue lines) and Y s (red curves) as a function of x = µ/T , for m s = 130 GeV, and m ψ = 120 GeV (left plot) and 110 GeV (right plot). As Fig. 5(left) suggests, Y s depends strongly on the t-channel semi-annihilation, as it can be compared the dashed  and solid red curves, with the former not considering the process in the Boltzmann equation and the latter containing it. As the mass difference between the two singlets increases, the effects becomes sharper, as shown in Fig. 5(right). Y ψ is almost independent on this process, as it can be seen through the overlapping of the dashed and solid blue lines. This effect is particularly interesting due to the fact that as Ω s becomes negligible, no direct detection bounds will apply on the pseudoscalar DM. This behavior was also seen in a two-component DM model consisting of two complex scalars stabilized by a Z 5 symmetry [54].
For ∆ > 0, the two semi-annihilations enter into the coupled Boltzmann equations, and it is not longer possible to assume Y ψ constant, therefore Eq. (12) is not a good approximation to determine Y ψ . In Fig. 6 we show the dependence of the relic abundance as a function of m ψ (left) and m s (right). Contrary to the case in the low mass regime of Fig. 4, in the left plot of Fig. 6 we observe that conversions occur at higher values of m ψ than s-channel semi-annihilations. When the latter opens up, the total relic abundance decreases in various order of magnitude. There is an interesting effect we want to point out, namely the fall of Ω s as m s m ψ . In both plots of Fig. 6 is possible to observe this effect, but in the right plot is more clear the role that conversions are playing, with the red dashed line making a tiny well from the soft growing as m s increases. This is, conversions of the type ss → ψψ start to be effective making a deviation from the well known behavior of the SHP [5] (for works with related behaviors see [36,55,56]). This type of wells, along with the Higgs resonance at m s ≈ m h /2, will be of particular importance to evade direct detection constraints.
In conclusion, we have analyzed the effectiveness of annihilations, conversions and semi-annihilations

Phenomenology
In this section we present the relevant experimental constraints on the model along with a full scan of the parameter space. Once obtained the allowed parameter space, we explore indirect signals and we set upper limits from it.

Experimental Constraints
When m s ≤ m h /2, the Higgs boson can decay into two s DM particles, with a decay width given by This contributes to the invisible branching rate Br(h → inv) = Γ inv /(Γ SM + Γ inv ), where the total decay width of the Higgs into SM particles is given by Γ SM = 4.07 MeV [57]. Experimental searches put stringent constraints on this quantity, and the most strict value is given by Br(h → inv) < 0.19 at 95% C.L. [57]. From the DM sector we have the following constraints. First, the measurement of the DM relic density today, given by Ω c h 2 = 0.120 ± 0.001 [58]. Based on the uncertainties in our computation, we apply this constraint with a tolerance of ∼ 10%, i.e. Ω c h 2 ∈ [0.11, 0.12]. Secondly, direct detection, which in multicomponent DM scenarios the interaction of DM with nucleus matter through the spinindependent (SI) cross section comes from weighting it by a factor which takes into account the relative abundance Ω i over the Planck measured Ω c , i.e., where σ SI,i is the spin-independent cross section. At tree level, it is only the pseudoscalar DM which interacts with nuclei via the Higgs portal, with its SI cross section given by [5] with m n as the nucleon mass, f N is a factor proportional to the nucleon matrix elements, and µ = m n m s /(m n +m s ) is the DM-nucleon reduced mass. Upper bounds on σ SI are given by XENON1T data [59] and the projections of XENONnT [60].

Scan
In this section we show a scan on the parameter space in the range m ψ,s ∈ [10, 1000] GeV (subject to m s < 2m ψ ) and g ψ , λ hs ∈ [0.001, 4π]. The full scan highlight interesting features that were already anticipated in Sec. 3.2. In Fig. 7 we show the scan of points projected on the plane (m ψ , m s ), with the density color representing the total abundance Ωh 2 (left), Ω ψ h 2 (middle) and Ω s h 2 , respectively. The green region (top left) corresponds to those points with m s > 2m ψ , then making s unstable. We pay attention to four regions based on Fig. 7 (left): (a) Dark red band. In general, Ω ψ get too big values for m ψ < m s , since no effective annihilation channels for the singlet fermion are present (see Sec. 3.2.1). However, there is a smooth relic density transition from the deep red to lower relic densities, due to the thermal tail distribution of conversions and s-channel semi-annihilation as ∆ 0. With respect to Ω s , it is always subabundant, specially for m s > m h where the t-channel semi-annihilations are present, and as we have described in Sec. 3.2.2, tiny mass shift between ψ and s makes Ω s to decrease strongly (deep blue region in Fig. 7 (right)).
(b) Green area in below. In this region the relic is mostly dominated by the pseudoscalar, as it can be seen from the plots. As it was shown in Fig. 4, Ω ψ < Ω s for conversion-driven processes of the type ψψ → ss and moderate couplings (g ψ ∼ 1). For very high m ψ is possible to observe that Ω ψ tend to increase, as conversions and semi-annihilations are less effective.
(c) Higgs resonance and thresholds. In the region 50 GeV m s ≤ m h the Higgs resonance and SM thresholds are present, then pseudoscalar annihilations are enhanced. These effects are more clear in Fig. 7 (right) for Ω s . For m s m h /2, s-channel annihilations of the singlet pseudoscalar via Higgs resonance tend to be very effective, decreasing Ω s substantially, as it can be seen with the deep blue horizontal line. For m h /2 < m s < m W the abundance of s lift up, and for m s m W the relic decreases again. The thresholds at m Z and m h are also present but they are less notorious.
(d) Blue-green upper region: For the singlet masses 100 GeV, the relic of both components tends to be low (blue points), but as the masses increase the effectiveness of conversions and semiannihilations become less strong, i.e. σv ∼ m −2 ψ,s (Appx. A), then raising the relic of both components. In this mass region the relic abundance seems to acquire moderate values in order to fulfill the Planck measurement.
The conjunction of the three plots suggests that most of the points that could give the correct relic abundance are just in very specific sectors, and those regions must to shrink even more when other constraints are taken into account. In Fig. 8 (left) we show the resulting parameter space points after imposing the constraints from relic abundance, Higgs to invisible and direct detection bounds. Specifically, the regions are: (i) Higgs resonance, (ii) region with m s > m h and ∆ < 0, and (iii) in the high mass with ∆ > 0 and m ψ < m s , indicating the corresponding Ω ψ with a color bar. Region (i) is the only one that allows m ψ > m s and where the Higgs to invisible constraint applies 3 . The points in region (ii) easily evade XENON1T constraints due to the fact that Ω s → 0, then the condition Ω ψ = Ω c is enough to fulfill all the experimental constraints. Finally, in the high mass region we found (right) The same selected random points shown in the left plot, now projected on (m s ,σ SI,s ), where we distinguish points with ∆ > 0 (orange) and ∆ < 0 (red). The continuous blue line corresponds to the upper limit given by XENON1T (1 t·y) [59], whereas the green dashed line a projection for XENONnT (20 t·y) [60].
some points anticipated by the analysis around Fig. 6, where the power of conversions and s-channel semi-annihilations makes Ω s to decrease enough evading XENON1T. The cost of the latter implies values for g ψ near the perturbativity limit criteria (this point is discussed in the last section). In Fig. 8 (right) we project the scan points on the plane (m s ,σ SI,s ), contrasted with the limits given by XENON1T and the XENONnT projection [60]. Even when XENON experiments rule out most of the orange points (i.e. regions (i) and (iii), a fraction of red points (∆ < 0 and m s > m h ) easily evade the strongest direct detection constraints. It is worth to say that this inverted peak at the Higgs resonance is a typical characteristic of Higgs portals, hardly to be ruled out by experiments.
To end this subsection, we would like to point out about the limitations of our random scan, which was based on overlaying exclusion limits (mainly from relic abundance and direct detection). It is well known that this type of exclusions present some drawbacks, such as the blindness to the influence of the choice of some parameters (e.g. DM halo distribution or the Higgs mass pole) on the experimental limits, or the fact that the allowed regions do not present additional information about which points are more favored than others. Considering those limitations, statistical analysis based on combined likelihoods functions may give a more accurate and realistic way to cure those disadvantages, although beyond the scope of the present analysis [8,62]

Indirect detection
In the previous subsection we found that three regions of the parameter space fulfill invisible Higgs decay, relic abundance and direct detection constraints, with two of them being practically in the region ∆ > 0. In that case, s-channel semi-annihilations are expected to give sizable fluxes of particles today through its s-wave nature. The t-channel semi-annihilation also goes in the s-wave, but we have checked that in general its cross section is lower than the corresponding s-channel, and indirect signals with a pair of s in the initial state have shown to be not sizable outside the Higgs resonance [5]. Additionally, Ω ψ h 2 reaches ∼ 50% to 99% of the total DM budget in the high mass regime (see Fig. 8(left)), consequently the scaling factor (Ω ψ /Ω c ) 2 does not considerable suppress the flux produced by the s-channel semi-annihilation. In the following, we show the box-shaped differential spectra [63] for this channel with its subsequent decay h → γγ, and secondly, restrictions on the parameter space coming from bounds based on searches of gamma rays (Fermi-LAT), anti-protons (AMS-02), and projections from the Cherenkov Telescope Array (CTA) are presented.

Box-shape gamma ray
The differential flux of photons produced in fermionic DM annihilations and received at earth from a given solid angle in the sky ∆Ω with a detector of area A is given by where σ sh v is the corresponding average annihilation cross section times velocity of the process ψ +ψ → s + h, Br h is the branching ratio of the Higgs into two photons, and dN/dE γ the corresponding normalized spectra. The J factor is the integral of the squared DM density ρ DM along the line of sight J = l.o.s dsρ 2 DM . We consider as our main region of interest the galactic center, which features ∆Ω = 1.30 sr, ∆Ω JdΩ = 9.2 × 10 22 GeV 2 cm −5 , assuming a NFW profile normalized to a local DM density of 0.4 GeV/cm 3 [64]. To determine the normalized spectra of emitted photons, we first note that their energies, in the fermion DM collision frame 4 , are given by where θ corresponds to the angle sustained by one of the photons and the in-flight Higgs. The energy of the emitted scalar particles in the rest frame of the collision of ψ are given by 4 Today, fermion DM ψ moves non-relativistically, therefore colliding practically in the earth rest frame.
For a fixed E h , the energy of the photon received at earth depends only on θ, with a maximum (minimum) energy for θ = 0(π/2) 5 then displaying a box-shaped spectrum centered at E c ≡ (E(0) + E(π))/2 = E h /2 and with a width ∆E ≡ E(0) − E(π) = E 2 h − m 2 h . Therefore, the normalized spectra can be written as 6 .
where the factor multiplicative factor 2 accounts for the two emitted photons (then E γ any of the two photons in 17), and Θ are Heaviside functions. In Fig. 9 (left) we show the box-shaped spectra for two random points allowed by the relic density abundance and XENON1T, contrasted with the data given by Fermi-LAT (red dashed line), where we have taken the background fitting function in order to estimate the data signal: dΦ γ /dE = 2.4 × 10 −5 (E γ /GeV) −2.55 GeV −1 cm −2 s −1 sr −1 [64]. As it was anticipated above, the differential fluxes for the selected points are well below the red dashed line, in principle not showing any possible tension. The blue and orange dot curves show the spectra when the two DM components become degenerated, widening the boxes and therefore peaking at higher energy photons. Only points with masses (and high couplings) as small as m ψ ∼ 100 GeV may approach significantly to the red dashed line, as it is depicted by the grey dashed line in Fig. 9 (left). It is known that the power of box-shaped gamma rays goes in their potential deviation with respect to the power-law Fermi-LAT background, then the previous result is premature to discard a possible constraint on the parameter space. In order to do this, in the following we use the recent upper bounds on gamma rays produced in semi-annihilations based on Fermi-LAT data, along with anti-protons flux measurements.

Upper bounds
Upper limits on σ hs v have not been constructed yet in order to constraint the parameter space via indirect searches. However, they can be inferred from existing upper bounds on the average annihilation cross section times velocity for DM + DM → DM + h process based on gamma-ray signals [66], and from DM + DM → bb(W + W − ) based on anti-proton flux measurements [67]. In the following we derive the useful algebraic relations to translate the existing upper bounds to our average semi-annihilation cross section times velocity σ sh v .
First, let us consider the process ψψ → ψ h, with ψ and ψ arbitrary DM particles, and h the : Differential gamma ray flux originated from the process ψ +ψ → s + h, followed by the decay h → γγ, for two random points fulfilling relic abundance and direct detection (continuous lines), whereas the dotted curves represent hypothetical points in which the two DM components are completely degenerated. The dashed red line corresponds to the flux given by Fermi-LAT [63], and the gray an hypothetical point with low (m ψ , m s ).
Higgs boson. This process presents an average cross section times relative velocity given by with J an arbitrary J-factor, m ψ the mass of the initial states, Φ h is the flux, and the outgoing Higgs having an energy E h (equivalently to 18). Now, it is possible to have the same flux Φ h with an outgoing Higgs having the same energy E h in another process given by ψψ → hh: where the 1/2 factor comes from the fact that we are considering one outgoing Higgs, and we have set m ψ = E h . Combining 20 and 21 we obtain with ξ x ≡ m 2 , with x being some particle in the final state. From the last relation ψ is an arbitrary state, then in eq. 22 taking ψ = ψ and ψ = s and combining the resulting expressions, Figure 10: Points with ∆ > 0 fulfilling direct detection constraints as a function of the re-scaled average annihilation cross section times velocity given by expressions 23 and 26 multiplied by the relative abundance (Ω ψ /Ω c ) 2 of the annihilating fermion DM. Upper bounds for Fermi-LAT [66], CTA [66] and AMS-02 [67] are given by the red continuous lines. For details see the main text.
it follows that Upper limits on σ ψh v [66], now can be compared to our σ sh v for certain values of m ψ and m s through eq. 23. In Fig. 10 (top row ) we show the projection of under-abundant points (green) and those that give approximately the correct relic abundance (blue points) for the resulting average annihilation cross section times velocity (left side of eq. 23) times the relative abundance of the annihilating DM particles. The red line in both plots corresponds to the upper bound found in [66] (right side of eq. 23). All these points are in the region ∆ > 0 and fulfill direct detection constraints.
None of the bounds touch the parameter space points, then not showing any possible exclusion. Note that the blue points well below the upper bounds (red line) correspond to the Higgs resonance region, since lower couplings (g ψ , λ hs ) are necessary to compensate the resonance effect, then reducing the ID signals.
On the other hand, AMS-02 anti-protons measurement set upper bounds on σ DM DM→XX v ≡ σ XX v with XX = bb, W + W − [67]. From the semi-annihilation with the Higgs decaying into bb we have that Additionally, the flux Φ bb in 24 can be produced by an arbitrary interaction ψψ → bb: Combining 24 and 25 we obtain Equivalently for W + W − . In Fig. 9 (bottom row ) we show the projection of points considering the semi-annihilation with Higgs decay into bb and W + W − , contrasted with the upper bounds given by AMS-02. From the resulting plots, bb search set the stringent bounds on the parameter space of the model, discarding points with the correct relic abundance with masses below ∼ 500 − 700 GeV, approximately.

Discussion and Conclusions
In this paper we have explored a simple extension to the SM containing two gauge-singlet fields: a Dirac fermion and a real pseudoscalar. From the DM point of view, the model present multiple scenarios with one or two DM components. Previous works have focused on the scenario when the Higgs portal coupling takes very small values, in such a way that the dark sector evolves decoupled from the SM, with the DM relic abundance produced via dark freeze-out/freeze-in. In the present work we found that as the Higgs portal take sizable values it is possible to produce a one-component DM candidate via freeze-in or a two-component DM via freeze-out. We have focused on the latter scenario, which only presents four free parameters: two masses and two couplings. The stability of the two singlets is guaranteed by parity arguments without the necessity of invoking additional symmetries. Furthermore, the relic abundance of both singlets is determined via freezeout through annihilations, DM conversions and semi-annihilations. The appearance of the latter is contrary to the standard belief that this type of processes appear only in the presence of symmetries larger than Z 2 . We have explored the relic density abundance of both DM components in the parameter space, founding interesting behaviors depending on the mass hierarchy and couplings values. Semiannihilations and DM conversions play an important role in two of the three available regions, making the pseudoscalar relic abundance low enough to evade direct detection XENON1T bounds, and then allowing DM with masses of hundreds of GeV up to the TeV scale. As it is usual with Higgs portals, the resonance region will not be discarded completely, even with the powerful projection of XENONnT.
We have complemented our analysis with indirect detection signatures and bounds. Firstly, we have explored the box-shape gamma rays signals appearing from a fermion DM semi-annihilation, showing small boxes signals for the fermion DM s-channel semi-annihilation. Furthermore, we have translated bounds from gamma-ray searches from Fermi-LAT and CTA projections onto the semiannihilation, not showing any possible tension, although the sensitivity of both experiments rely in the ballpark of part of the parameter space of the model. In fact, considering that our numerical methods are not very precise, a more detailed and sophisticated analysis such as a global-fit could be perfectly sensitive to these gamma-rays upper bounds in view of the expected quantitative differences between the two methods. As a final analysis, we have tested the model with anti-protons upper bounds from AMS-02 experiment, showing an exclusion of fermion DM masses below ∼ 500 − 700 GeV.
Finally, considering that the viability of the model requires high values for the dark sector coupling g ψ in some regions of the parameter space, there are two important points to be taken into account, although the precise calculation of them are beyond the scope of the present work. The first is related to the possible appearance of Landau poles at low energy scales, implying an urgent UV completion of the model (e.g. embedding the pseudoscalar into a complex singlet which acquires vev). The second point is related to the spin-independent one-loop amplitude in direct detection for the singlet fermion, which due to the high coupling values it may give a sizable contribution, possibly excluding even more parameter space of the model, in particular the special region in which the relic abundance of the pseudoscalar drops to zero.
with m 1 and m 2 the masses of the annihilating particles, keeping the s-wave only when it is present: The corresponding expressions for the annihilations of a pair of s into SM particles can be found in the appendix of [5].