Bouncing pNGB dark matter via a fermion dark matter

In addition to the Standard Model, the introduction of a singlet complex scalar field that acquires vacuum expectation value may give rise to a cosmologically stable pseudo-Nambu-Goldstone boson (pNGB), a suitable dark matter (DM) candidate. This work extends this scenario by including a second cosmologically stable particle: a fermion singlet. The pNGB and the new fermion can be regarded as DM candidates simultaneously, both interacting with the Standard Model through Higgs portals via two non-degenerate Higgs bosons. We explore the thermal freeze-out of this scenario, with particular emphasis on the increasing yield of the pNGB before it completely decouples (recently called Bouncing DM). We test the model under collider bounds, relic abundance, and direct detection, and we explore some indirect detection observables today.


Introduction
The dark matter (DM) field is an active area of research comprehending new extensions to the Standard Model (SM) of particle physics.One simple well-motivated scenario is the global U (1) complex scalar extension to the SM, which after the spontaneous and explicit symmetry breaking gives rise to a pseudo-Nambu-Goldstone boson (pNGB) and a second Higgs boson [1,2,3].This pNGB, being cosmologically stable, has been studied as a thermal relic in a standard freeze-out, presenting viable DM masses in the ballpark of the electroweak scale.A second Higgs may be relevant in the explanations of anomalies seen at LEP and LHC [4].
On the other hand, the question of whether the DM genesis was thermal or not has given rise to interesting possibilities.For instance, non-thermal DM, such as feebly-interacting massive particles (FIMPS) [5], or DM subject to an exponential growth [6], evolve increasing their yields from a negligible initial abundance until Hubble expansion dominates and the relic abundance is set.Analogously, could thermal relics present a feature like this growing yield in the early universe?Recently, it has been suggested that in specific scenarios, the yield of a DM particle could have been undergoing a different pattern of decoupling, with a period of exponential growth before complete freeze-out.This idea has been called bouncing dark matter [7].Although this bouncing effect in the yield of some species has been already observed in the context of thermal DM, see e.g.[8,9,10,11], the connection of this feature with indirect detection has only been worked out in [7].
In this work, we study a multi-component DM scenario in which one of the stable relics undergoes this bouncing effect.The model consists of the introduction of two gauge singlets: a fermion and a complex scalar field, both transforming under a global symmetry 1 .After spontaneous symmetry breaking (SSB), the singlet scalar gets vacuum expectation value, giving rise to a second Higgs boson which will mix with the Higgs excitation of the SM Higgs doublet, and a Nambu-Goldstone boson χ.Imposing a U (1) soft-breaking term in the scalar sector, χ becomes a cosmological stable pNGB if its mass is lower than twice the fermion singlet mass.In this way, the model presents two DM components that communicate to the SM through two Higgs-like particles.
We study the thermal freeze-out of both stable relics, with special emphasis on the details of the bouncing effect undergone by the yield of the pNGB, and also some prospects of indirect detection.Finally, as the minimal scenario of the pNGB as the only source of DM allows masses above 50 GeV [3], unless resonance annihilation effects are present, in the present multi-component scenario we explore the viability to have pNGB DM masses below 50 GeV.
The paper is organized in the following way.In Sec. 2 we introduce the model.In Sec 3 we present the Boltzmann equations for the system, with a detailed analysis of the relic abundance of both components, highlighting the bouncing of the pNGB yield, and some average cross sections relevant for indirect detection.In Sec. 4 we present the phenomenology of the model for two-DM components, imposing relevant constraints, and studying some indirect signals.Finally, in Sec 5 we discuss and conclude our work.

Model
Aside from the SM particle content, we add two SM singlets, a Dirac fermion ψ (for the Majorana case the construction follows equivalently, e.g.[13]) and a complex scalar S. Decomposing ψ = ψ L + ψ R , we impose a chiral approximate global symmetry U (1) V × U (1) A , with the new fields transforming as [16] U with β V,A arbitrary constants.More concisely, the spinor under U (1) V × U (1) A transform as ψ → e i(β V +γ 5 β A ) ψ.This symmetry and particle content give rise to the Lagrangian: with the potential given by Here H is the Higgs field, and the U (1) A soft breaking term is given by where we have assumed that a subgroup Z 2 of U (1) A remains unbroken (i.e.β A = π), in such a way that only soft-breaking terms with even powers of S are allowed.The couplings g ψ and m 2 χ can be made real by field redefinitions.Note that the global symmetry does not allow a mass term for the fermion singlet, but it will be generated through the spontaneous symmetry breaking of S. Analogously to the global U (1) complex scalar model [3], the Lagrangian is invariant under a CP-symmetry: ψ L,R → ψ R,L and S → S * .Notice that Z 2 symmetry in which S → −S is broken by the vacuum v s [17].
The stationary point conditions at (h, s) = (0, 0) are As µ 2 S < 0, S triggers the SSB, making U (1) A symmetry is spontaneously broken.Explicitly, choosing the Hermitian basis, S = (s ′ + iχ)/ √ 2, the minimum is acquired when s ′ 2 + χ 2 = v 2 s > 0, with v s the vacuum expectation value of the complex field.Without loss of generality, we set ⟨s ′ ⟩ = v s and ⟨χ⟩ = 0. Noting that ψL ψ R + ψR ψ L = ψψ and ψL ψ R − ψR ψ L = ψγ 5 ψ, the interactions around the new vacuum after the SSB of U (1) A are given by Note that the fermion has acquired a mass m ψ ≡ g ψ v s / √ 2, a scalar interaction with s, and a pseudo-scalar interaction with χ.After electroweak symmetry breaking (EWSB), the Higgs field mixes with s, which introduces a mixing angle θ which satisfies with v h = 246 GeV, h = cos θh 1 + sin θh 2 and s = − sin θh 1 + cos θh 2 , and the masses of the scalars are where we identify h 1 with the 125 GeV Higgs boson.From this last relation, we see that In the physical basis, the interaction of the singlet fermion with the scalars is given by with the potential V (h 1 , h 2 , χ) explicitly detailed in the App.C. The free parameters of the model are three masses and two couplings: The model described by eq. 12 may present one or two DM candidates, depending on the mass hierarchy between ψ and χ.If m χ ≥ 2m ψ , the fermion is the only stable field, and for m χ < 2m ψ the pNGB becomes stable, then the model presents two DM candidates (for a more detailed study of the stability of the pNGB, see [18].Also see [19]).In this work, we consider the latter case.
Figure 1: Tree level diagrams relevant for the calculation of the relic abundance of ψ and χ.Here X refers to an SM particle, including h i with i = 1, 2.

Relic abundances and cross sections
In the following, we analyze the relic abundance of the two-component DM scenario previously described obtained via freeze-out of each DM species, analyzing the hierarchy between the DM particles in order to distinguish the bouncing effect.Furthermore, we analyze cross-sections at low temperature relevant for indirect detection observables.

Boltzmann equations
We assume that in the early Universe both DM candidates were in thermal equilibrium with the SM particles.In Fig. 1 we show the corresponding Feynman diagrams that participate in the relic density calculation.Based on [20] and without loss of generality assuming that m ψ < m χ , the evolution of the individual singlet abundances Y i ≡ n i /s, with i = ψ, χ, as a function of the temperature x ≡ µ/T , with µ = m ψ m χ /(m ψ + m χ ), are given by where we have defined where X here represents an SM particle.The entropy density s and Hubble rate H in a radiationdominated universe are 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 , respectively.The equilibrium densities, Y i,e ≡ n ie /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 equation gets modifications as the hierarchy of the singlets changes, and in this way, the equation may be derived using the detailed balance principle n a,e n b,e ⟨σ abcd v⟩ = n c,e n d,e ⟨σ cdab v⟩, with the indexes a, b, c, d for the respective particles.We implemented the model into LanHEP [21] and into micrOMEGAs code [22] to perform the calculations.For the rest of the paper, sometimes it will be convenient to use the quantity ∆ i = 2m ψ − m χ − m hi .

Mass hierarchies
We distinguish two hierarchies between the DM particles which can make the yield behavior quite different, particularly for χ: (i) m ψ > m χ , what we call the normal hierarchy, and (ii) m ψ < m χ , the inverse hierarchy.In the former case, the freeze-out of the DM particles results to follow the standard freeze-out of two interacting DM particles, each one decoupling from the SM plasma at x ≈ 15 − 20.We have exemplified this with a few parameter space points in Fig. 2(left) where we show the yields of each DM candidate as a function of the temperature.As the two DM particles present several interactions, they continue decreasing their yields for some time after each particle completely freezes out.
On the other hand, in the inverse hierarchy, and at high temperatures, ψ and χ being in thermal equilibrium with the SM plasma, we assume that both particles have vanishing chemical potentials, µ ψ = µ χ = 02 .Since m χ > m ψ , the fermion is less Boltzmann suppressed than χ, and assuming that ψ does not rise in abundance and they keep the same chemical potential, one has that n χ = (n χ,e /n ψ,e )n ψ ∼ e −(mχ−m ψ )/T n ψ , i.e., the abundance of χ decreases as T decreases.After chemical decoupling from the SM sector, the DM particles may develop a chemical potential n i ≈ n i,e e µi/T , in such a way that the yield of χ now results in  Provided µ χ > µ ψ and even having a decreasing n ψ , the Boltzmann suppression in eq.19 from the factor n χ,e /n ψ,e can be compensated by the exponential e (µχ−µ ψ )/T , rising exponentially the number density of χ for some time until Hubble expansion dominates.After the dark particles decouple chemically from the SM thermal bath, the process sustaining chemical equilibrium within the dark sector are semi-annihilations ψ ψ ↔ χh i , with i = 1, 2, such that µ χ ≈ 2µ ψ .The resulting effect is shown in Fig. 2 (middle) for a specific choice of parameters, where the rising of Y χ results for a finite interval of temperature.The size of the yield increment is model-dependent, and as it can be noted in Fig. 2 (middle), the height of the bouncing depends on m h2 .The decreasing in the yield of the bouncing comes from the kinematical suppression of the process responsible for it, i.e. ψ ψ → χh 2 , since 2m ψ < m χ + m h2 , suppressing ⟨σ ψψχh2 v⟩, in turn decoupling χ and ψ at earlier times in comparison to the case with lower m h2 .We say that for those cases we have ∆ 2 < 0 (for this particular case, we also have ∆ 1 < 0)3 .In Fig. 2 (right), we observe graphically the non-zero values taken by the chemical potential of the two stable particles, fulfilling µ χ ≈ 2µ ψ , with the highest value of µ χ − µ ψ being the case with the strongest bouncing.To obtain each chemical potential, we have solved µ i from n i ≈ n i,e e µi/T , such that with Y i (T ) obtained with micrOMEGAs code after solving the coupled Boltzmann equations (cBE), eq.14.
Another relevant feature regarding the bouncing is the variation in the minimum Y χ at which the bouncing starts to appear.If the pNGB remains more time in thermal equilibrium with the SM, the bouncing will start at later times, and it may grow the yield several orders of magnitude before coming to an end.This is exemplified in Fig. 3(left) for the parameters indicated in the plot.As shown, in this case, it is clear that the bigger m χ , the bigger the bouncing.This effect is due to the fact that the semi-annihilation term is proportional to λ ψ ψχhi and impacts directly  in the decoupling temperature of χ from the thermal bath, with a strong dependence on m χ .In App.A we study this effect in greater detail from a numerical and semi-analytical way.
As explained before, the bouncing effect proceeds through either semi-annihilation ψ ψ → χh 1 or ψ ψ → χh 2 , each process with a dependence of sin 2 θ and cos 2 θ, respectively, affecting the relic abundance directly.For instance, in Fig. 3(right) we plot the relic abundance as a function of the second Higgs mass, for m ψ = 1000 GeV, m χ = 1500 GeV and g ψ = 1, for different values of θ.For this parameter choice, the bouncing effect is present for all the shown combinations of values of (m h2 , θ).As it is shown by the solid red line (θ = 0.1), Ω χ remains constant for m h2 ≳ 500 GeV since, even though ⟨σv⟩ ψ ψχh2 is kinematically suppressed (i.e., the pink region representing ∆ 2 < 0), the process ψ ψ → χh 1 becomes the leading one for the bouncing effect.On the contrary, for a more decoupled dark sector-SM case, e.g.θ = 10 −5 , ψ ψ → χh 1 process becomes suppressed by sin θ, and ψ ψ → χh 2 remains as the effective process, with Ω χ decreasing as m h2 gets higher values 4 .In that case, the bouncing effect becomes smaller.Therefore, the mixing angle and m h2 are important parameters behind the bouncing effect, impacting the value of Ω χ , and also determining the leading semi-annihilation in the relic density calculation.

Cross sections
In this subsection, we study the values of some of the most relevant average annihilation cross section times relative velocity at v rel → 0 from DM (semi)annihilation, especially when they surpass the thermal canonical value 2 − 3 × 10 −26 cm 3 /s [7].We run two random scans, one considering m χ = 300 GeV and the other m χ = 700 GeV.In both scans we have considered m ψ = 500 GeV, m h2 = [50, 2000] GeV, g ψ = [0.1,10] and tan θ = [10 −2 , 10 1 ].We have selected all the points which match the observed relic abundance.As shown in the first row of Fig. 4, we have projected the points in different planes, with the color of each point indicating the corresponding value of tan θ.As it is expected in the normal hierarchy m ψ > m χ , the first two plots corroborate the fact that the weighted cross sections never surpass the thermal cross section reference (pink regions) 5 .The other two plots in the right of the first line show the corresponding values for the couplings of the model.
In the lower row of Fig. 4, we project the random scan for the inverse hierarchy m ψ < m χ .Contrary to the previous case, here many points surpass the thermal reference value in the first two plots, although with clear differences in their value acquired by θ.Notice that for the case of fermion DM annihilation, those points above the thermal relic tend to be disfavored by the high values of g ψ , due to perturbativity.For the case of the pNGB DM annihilation, the points with strong cross-sections seem to respect perturbativity for g ψ and λ HS .Nonetheless, as we show in the next section, θ ≳ 0.1 enter in conflict with both direct detection and collider constraints.
Finally, a few comments to highlight.First, analogously to the 125 GeV Higgs boson h 1 , we have been assuming that the mediator h 2 remains in thermal equilibrium with the SM plasma at any moment, so its chemical potential vanishes.Secondly, a necessary condition for the existence of the bouncing is m χ > m ψ , but it is also present in either case m ψ > m h2 or m ψ < m h2 .Third, ∆ 1 > 0 and/or ∆ 2 > 0 do not guarantee the bouncing effect.In the following section we study the phenomenology of the model.

Phenomenology
In this section we consider the setup at hand as a fully realistic DM model, considering constraints from collider and DM searches.We explore the viability of pNGB DM below 50 GeV, and also we study the normal and inverse hierarchy for masses of hundreds of GeV.Finally, we explore a few indirect detection observables.

Dark matter constraints
The relic abundance measure today is given by the most updated Planck result [27].We take the value Ω c h 2 = 0.12 as the measured DM relic density, with an error of 10% in the calculation of it with micrOMEGAs.As usual, the total relic abundance in the two-component DM model is given by the sum of each stable particle, i.e.Ωh 2 = Ω ψ h 2 + Ω χ h 2 .
On the other hand, direct detection (DD) for the pNGB DM candidate is relaxed due to its Goldstone nature [3] (unless new explicit global symmetry-breaking sources are present [28]).It has been shown that the direct detection rate of a pNGB DM particle is too low to be observed with present experiments [3], even at the one-loop level [29], then we do not take into account the relative contribution of the pNGB for DD constraints 6 .On the contrary, ψ is subject to sizable constraints appearing from the scattering in t-channel exchange of h 1 and h 2 .At tree level, it is given by [13] where m N denotes the nucleon mass and f p = f n ≈ 0.27.Even though in most of the cases we take the most recent bounds on direct detection from Lux-Zeplin experiment (LZ) [30], in some cases we also used XENON1T [31] and XENONnT projection [32].As an example of the intensity of LZ constraints on the parameter space of the model, in Fig. 5(left) we show the resulting allowed parameter space for a representative fermion mass of 500 GeV, taking g ψ = π (black lines) and g ψ = 0.5 (grey lines).The region above each curve is excluded by LZ, and the corresponding dashed lines represent the bound when the fermion acquires 10% of the total relic abundance.As it can be appreciated, direct detection bounds are strong, particularly for m h2 < m h1 .

Collider constraints
A second Higgs is constrained throughout the combination of its mass and coupling to the SM particles, i.e. its mixing angle.Direct searches of a second Higgs set that θ ≲ 0.15 for m h2 < 50 GeV, whereas for m h2 > 100 GeV, electroweak precision tests (EWPT) set that θ ≲ 0.3 [33,34].
If either of the DM particles has a mass below half of the mass of the 125 GeV Higgs boson, the latter may decay into a pair of DM particles, resulting in an invisible decay of the Higgs boson at colliders.Measurements set limits on the branching fraction of the Higgs into invisible particles, with the most updated value being BR(h 1 → invisible) ≲ 14% at 95% C.L. [35].In our scenario, if the kinematic is allowed, we have that where each decay width is given by with r ≡ m h1 /m ψ .For the parameter space that we are interested in, h 2 is short-lived7 , therefore the decay h 1 → h 2 h 2 gives not invisible products, but it is important for the total width decay of the Higgs boson.Furthermore, as h 2 is short-lived, it does not affect constraints from Big Bang Nucleosynthesis (BBN).Finally, electroweak precision test set constraints on the combination (m h2 , θ), but for m h2 < m h1 /2 they result to be much weaker than Higgs to invisible constraints [33].

pNGB mass below m h /2
In the original scenario of pNGB DM [3], masses for the latter below ∼ 50 GeV are normally not possible, since Higgs to invisible bounds on the Higgs portal coupling become relevant, and to fulfill the correct relic abundance, the Higgs portal can not take arbitrarily small values, otherwise an overabundance is obtained.One way to avoid the latter is having 2m χ ∼ m h2 in such a way that the annihilation of the pNGB occurs on-shell in the s-channel, then decreasing sufficiently the Higgs portal coupling to evade the Higgs to invisible constraint (see e.g.[4]).In our twocomponent DM scenario, we show that it is possible to have viable pNGB with masses below m h /2 without entering into resonance effects, fulfilling the correct relic abundance, and evading Higgs to invisible bounds and direct detection.The dynamics for the calculation of the relic abundance is independent of the mixing angle provided θ ≲ 0.1, since all the diagrams containing h 1 are suppressed by tan 2 θ, then the relic abundance is determined by the fields of the dark sector and h 2 .In Fig. 5(middle) we show the relic abundance behavior for the two DM candidates as a function of m ψ , assuming m χ = 30 GeV, m h2 = 29 GeV, g ψ = 1 and θ < 0.1.Here, m ψ ≈ 55 GeV is the required value to achieve the correct relic abundance.Notice that m h2 must be lighter than the two DM components, in order to avoid overabundance.In the same line, we observe in the right plot of Fig. 5, that the maximum value that θ can take to evade LZ is ∼ 0.006.Notice here that the solid line considers the scaling of σ SI by the fraction of relic abundance of the fermion.It is interesting to observe here that constraints from 125-GeV Higgs decaying into invisible particles are complementary to direct detection (green regions), although less strong than LZ in this case.Notice that future projections from XENONnT and collider searches will be competitive between them, excluding even much more parameter space.GeV and g ψ = 1 (i.e., the parameter space point that in the plot in the middle fulfill the correct relic abundance).In red we show constraints from XENON1T, LZ and the projection of XENONnT, and in green the excluded regions by the invisible Higgs decay considering a branching of 0.14 (dark green) and a projection of 0.01 (light green).In summary, it is possible to have pNGB DM for masses below ∼ 50 GeV, without recurring to resonance effects.The price to pay in order to evade strong direct detection and Higgs to invisible constraints is to decrease θ enough in such a way to compensate the growing behavior of the spin-independent cross section with light h 2 , since σ SI ∼ sin 2 (2θ)/m 4 h2 .

Indirect detection prospects
Considering that the two-DM scenario presents sizable average annihilation cross section at low temperatures, now we focus on specific (semi)annihilations prospects relevant for indirect detection observables.We pay special attention to the (semi)annihilation processes ψ ψ → χh i and χχ → XX, with X an SM state8 .The corresponding partial average annihilation cross-sections have been calculated using CalcHEP, expanding around v rel = 0, in order to retain the s-wave contribution only (see App. B for the resulting analytical expressions).
The partial cross sections result to be highly dependent on the parameters of the model.As we exemplify in the upper row of Fig. 6, considering m ψ = 500 GeV, m h2 = 600 GeV and g ψ getting the appropriate value to match the correct relic abundance, the cross sections not only vary by orders of magnitude depending on m χ , but as θ decreases, the parameter space available to obtain the correct relic abundance shrinks allowing only m χ ≈ m h2 /2, otherwise an overabundance is obtained.In this way, in the normal hierarchy and small mixing angles it is possible to obtain sizable cross sections, as shown by the orange solid line and the dashed curves in the left plot of the upper row of Fig. 6, but only in a reduced parameter space.On the contrary, higher tan θ values,  In each plot we have taken m ψ = 500 GeV, m h2 = 600 GeV, and tan θ = (10 −3 , 10 −2 , 10 −1 ) (from left to right).g ψ takes the necessary value to obtain the correct relic abundance, considering perturbativity constraints.The pink band represents the thermal canonical value 2 − 3 × 10 −26 cm 3 /s.The solid lines are the values corresponding to the fermion DM annihilation, whereas the dashed lines are the corresponding pNGB annihilation.(bottom) Zero-velocity average annihilation cross section for different DM (semi)annihilation channels as a function of the mixing angle.The values of the parameters here have been taken as m ψ = 500 GeV, m χ = 800 GeV, m h2 = (130, 300, 600) GeV (from left to right in the plots), and g ψ takes the necessary value to obtain the correct relic abundance.The pink band is a reference for the thermal value, the grey one represents LZ bounds, and the red one collider constraints.e.g.tan θ = 0.1, imply less suppression for (semi)annihilation processes into SM states including h 1 in the final state, presenting strong cross sections specially in the case m χ > m ψ , where the bouncing effect is present.This can be seen in the third plot of the first row of Fig. 6, showing a wider range of m χ allowed.For completeness, we also present the case with tan θ = 10 −2 .
We confront the resulting zero-velocity relic weighted cross sections with direct detection bounds from LZ for three scenarios in which we vary m h2 , with each case fulfilling the correct relic abundance.In Fig. 6(below), we show the results as a function of the singlet-doublet mixing angle assuming m ψ = 500 GeV, m χ = 800 GeV, and m h2 = (130, 300, 600) GeV (from left to right, respectively).We have taken small values for m h2 to see the relaxing effect on direct detection bounds.LZ data rule out the shaded region in each plot.In the left plot of Fig. 6(below), LZ bounds are relaxed, due to the algebraic cancellation between the close-in mass of h 1 and h 2 .As m h2 deviates away from m h1 , LZ bounds start to be notorious and strong, as it is shown by the middle and the plot in the right, even for θ ≪ 1.In this way, all the cross sections with values above the thermal value obtained in the case in which m ψ < m h2 < m χ , result disfavored by LZ.
Notice also that if we take m h2 ≪ m h1 , direct detection becomes even stronger, therefore we do not look into that region.In any case, even though direct detection constrains the average crosssection above the thermal value (strongest bouncing effect regime), the model remains viable, with indirect detection near the ballpark of the thermal reference value.A detailed analysis of the full viable parameter space is beyond our goal in this work.

Conclusions
We have studied a simple extension to the SM which under reasonable assumptions presents two DM candidates simultaneously, a fermion and a pNGB, with each one communicating to the SM via the Higgs portal through two Higgs-like bosons.Assuming a thermal scenario in which both DM candidates freeze-out in a radiation-dominated universe, we have found a peculiar yield behavior for the pNGB when this is heavier than the fermion singlet: it bounces.This exponential yield increment may reach several orders of magnitude.This model is one of the first scenarios in which the bouncing effect is exemplified in more detail.
This effect has a direct impact on indirect searches.In this way, we have explored the zero velocity average annihilation cross sections relevant for indirect searches, finding parameter space regions in which both the fermion semi-annihilation and the pNGB annihilation today present values above the canonical thermal value.We have seen that as the fermion is subject to strong tree-level spin-independent DD constraints, the model requires a suppression on the mixing angle θ unless the two Higgses become too degenerated.In either case, collider, and specially LZ upper bounds, force to reduce some of the parameters of the model in such a way that the strongest indirect detection signals must be suppressed when they fulfill the correct thermal value.
Last but not least, we have shown that it is possible to recover pNGB DM for masses below 50 GeV, in contrast to the simple model in which invisible Higgs decays ruled out the parameter space for low masses (unless a resonance effect is present).The cost of this is to decrease the singlet-doublet mixing angle to values much below the unity.All in all, the model, being a multicomponent DM scenario, is not only elegant in its construction, but it presents the interesting effect of pNGB yield bouncing in the early universe, although the model itself suffers severe constraints from direct detection experiments. of each equation, solving them numerically and in a semi-analytical way.As we are interested in the bouncing effect details, we always consider m χ > m ψ .For the sake of simplicity, we define λ A1 ≡ λ ψψXX , λ A2 ≡ λ χχXX and λ S ≡ i=1,2 λ ψ ψχhi , then the cBE becomes First, we solve the cBE eq. 25 and 26 numerically in Python with solver ivp, assuming thermal equilibrium at x = 1.In order to solve the system, we assume m ψ = 100 GeV, and m χ = 150 (solid red) and m χ = 180 GeV (dashed red).The results are shown in the plot in the left of Fig. 7, where the solution for Y χ results to be highly sensitive to m χ .Here, we are assuming λ A1 = λ A2 = λ S , where each average annihilation cross section is taken to the value 10 −9 GeV −2 .As it is clear, the heavier is χ, the longer it stays in thermal equilibrium with the SM.This result has also been checked with micrOMEGAs.For comparison, in Fig. 7(middle plot), we show the resulting behavior considering λ S = 0, a typical freeze-out of two non-interacting DM particles.Thereby, the presence of the semi-annihilations in the cBE, results in a significant impact on Y χ , keeping χ longer in thermal equilibrium, with a strong dependence on m χ .
In the following, we take a semi-analytical approach to solve the cBE.Based on what we have gotten in the numerical solution, we assume that Y ψ tracks its equilibrium yield for x < 15, which is our temperature region of interest, then we take Y ψ = Y ψ,e .Furthermore, following the freeze-out approximation [37], we take Y χ ≈ (1 + δ)Y χ,e , with δ a positive number that grows slowly.After some algebra, and taking dδ/dx ≪ δ, we obtain where we have evaluated all the quantities at x f .This is a transcendental equation for x f which can be solved easily.We solve this equation considering δ f = 1 9 , which is the moment at which starts the chemical decoupling.As it can be seen in Fig. 7(right), the orange line represents the left side of eq.27, whereas the rest of the lines correspond to the r.h.s. of eq.27, with each line considering (w/s) and not (n/s) the semi-annihilation term proportional to λ S .As it was anticipated by the numerical solution in the first part of this appendix, the overall effect of the presence of semi-annihilations in the cBE, makes χ to be longer in thermal contact with the plasma, since the contribution of the second term in the r.h.s of eq.27 add a new positive contribution.Secondly, higher values of m χ result in a later chemical freeze-out temperature in the case in which semi-annihilations are present (w/s case).Therefore, the numerical and the semi-analytical approach for the temperature decoupling of Y χ agree, with the semi-annihilation not only predicting a bouncing as m χ > m ψ , but impacting strongly the behavior of Y χ before the bouncing.GeV −2 .The plot in the left considers the semi-annihilation term proportional to λ S , whereas the plot in middle does not.(right plot) Solution to the transcendental equation eq.27, with the intersection of the orange curve with the rest at x f .In the legend, each pair of values in the parenthesis represents (m ψ , m χ ) GeV.

B Cross sections
In this appendix, we present algebraic expressions for each DM average annihilation cross section times relative velocity relevant for the calculation of indirect detection.All the cross sections here were obtained with CalcHEP.As we expand each cross-section in powers of v n , with v the relative velocity of the colliding non-relativistic DM particles, with n taking positive even numbers, these expressions are not precise enough around poles or thresholds.In this way, we take the nonrelativistic expansion s = 4m 2 i 1 + v 2 4 , with i = ψ or χ, retaining the s and p wave only.For the case of the fermion DM, we have The expression for ⟨σv⟩ ψ ψχh2 is the same as the previous one but without tan 2 θ in the numerator, and m h1 → m h2 .Additionally, the annihilation of the fermion DM into gauge bosons is velocity suppressed: The annihilation cross-section for the pNGB are The expressions for ⟨σv⟩ χχhihj , with i, j = 1, 2, are too long to be written here.However, they result to be s-wave, therefore relevant for indirect detection observables.

C Scalar potential
The complete scalar potential takes the form:

Figure 2 :
Figure 2:In the first two plots we show the yields of ψ (solid lines) and χ (dashed lines), assuming g ψ = 1 and θ = 0.1.In the plot in the left we consider the normal hierarchy with (m ψ , m χ ) = (200, 150) GeV, whereas in the plot in the middle, the inverse hierarchy with (m ψ , m χ ) = (100, 150) GeV.The plot on the right shows the evolution of the chemical potentials of each DM particle.

Figure 4 :
Figure 4: Random scans with each point fulfilling the correct relic abundance, assuming m ψ = 500 GeV, and the color indicating the value of θ.In the first row we consider m χ = 300 GeV, whereas in the second one m χ = 700 GeV.In the first two y-axes, f i ≡ Ω i /Ω c , with Ω i the relative abundance of the i = ψ, χ initial state.For more details of the scan see the text.

Figure 5 :
Figure5: (left) Direct detection constraints on the fermion DM considering LZ bounds.The solid black (grey) curves represent the contours of σ SI for m ψ = 500 GeV and g ψ = π(0.5),whereas the respective dashed curves consider that the fermion contributes with a 10% of the total relic abundance.LZ data rule out the region above each curve.(Middle) Relic abundance in the lowmass regime, assuming m χ = 30 GeV, m h2 = 29 GeV, g ψ = 1, and θ < 0.1.The horizontal pink region corresponds to the observed relic abundance.(right) Direct detection and invisible Higgs decays constraint assuming m ψ = 55 GeV, m χ = 30 GeV, m h2 = 29 GeV and g ψ = 1 (i.e., the parameter space point that in the plot in the middle fulfill the correct relic abundance).In red we show constraints from XENON1T, LZ and the projection of XENONnT, and in green the excluded regions by the invisible Higgs decay considering a branching of 0.14 (dark green) and a projection of 0.01 (light green).

Figure 6 :
Figure 6: (top) Relic weighted average cross section times relative velocity as a function of m χ .In each plot we have taken m ψ = 500 GeV, m h2 = 600 GeV, and tan θ = (10 −3 , 10 −2 , 10 −1 ) (from left to right).g ψ takes the necessary value to obtain the correct relic abundance, considering perturbativity constraints.The pink band represents the thermal canonical value 2 − 3 × 10 −26 cm 3 /s.The solid lines are the values corresponding to the fermion DM annihilation, whereas the dashed lines are the corresponding pNGB annihilation.(bottom) Zero-velocity average annihilation cross section for different DM (semi)annihilation channels as a function of the mixing angle.The values of the parameters here have been taken as m ψ = 500 GeV, m χ = 800 GeV, m h2 = (130, 300, 600) GeV (from left to right in the plots), and g ψ takes the necessary value to obtain the correct relic abundance.The pink band is a reference for the thermal value, the grey one represents LZ bounds, and the red one collider constraints.

Figure 7 :
Figure 7: (left and middle) Numerical yield behavior for the cBE in eq. 25 and 26, assuming m ψ = 100 GeV, λ A1 = λ A2 = λ S , with each cross section set to the canonical value ⟨σv⟩ = 10 −9GeV −2 .The plot in the left considers the semi-annihilation term proportional to λ S , whereas the plot in middle does not.(right plot) Solution to the transcendental equation eq.27, with the intersection of the orange curve with the rest at x f .In the legend, each pair of values in the parenthesis represents (m ψ , m χ ) GeV.