A Minimal Model For Two-Component FIMP Dark Matter: A Basic Search
S. Peyman Zakeri S. Mohammad Moosavi Nejad Mohammadreza Zakeri S. Yaser Ayazi
hep-ph

	'α'
	0.65 cm
	

		 December 30, 2023



A Minimal Model For Two-Component FIMP Dark Matter: A Basic Search  0.5 cm S.  Peyman  Zakeri^1, S. Mohammad Moosavi Nejad^1,2, Mohammadreza  Zakeri^3,4 and S. Yaser Ayazi^50.5 cm
	
	^1Faculty of Physics,  Yazd University, P.O. Box 89195-741, Yazd, Iran0.1 cm^2School  of  Particles  and  Accelerators,  Institute  for  Research  in 
		Fundamental Sciences (IPM), P.O. Box  19395-5531, Tehran, Iran
		^3Physics and Astronomy Department, University of California, Riverside, California 92521, USA
		^4CAS Key Laboratory of Theoretical Physics, Institute of Theoretical Physics, Chinese Academy of Sciences, Zhong Guan Cun East Rd.55, Beijing, 100190, China
		^5Physics Department, 		Semnan University, P.O. Box 35131-19111, Semnan, Iran
	
	





	
In the multi-component configurations of dark matter phenomenology, we propose a minimal two-component configuration which is an extension of the Standard Model with only three new fields; one scalar and one fermion interact with the thermal soup through Higgs portal, mediated by  the other scalar in such a way that the stabilities of dark matter candidates are made simultaneously by an explicit Z_2 symmetry. Against the most common freeze-out framework, we look for dark matter particle signatures in the freeze-in scenario by evaluating the relic density and detection signals. A simple distinguishing feature of the model is the lack of dark matter conversion, so the dark matter components act individually and the model can be adapted entirely to both singlet scalar and singlet fermionic models, separately. We find dark matter self-interaction as the most promising approach to probe such feeble models. Although the scalar component satisfies this constraint, the fermionic one refuses it even in the resonant region.  
      








§ INTRODUCTION
 


Weakly interacting massive particles (WIMPs) are the most popular solution to the puzzle of dark matter (DM) <cit.>. In TeV scale (LHC scale) new physics, DM particles follow the thermal scenario in which they reach thermal and chemical equilibrium with the bath particles but lose it at the freeze-out temperature (which is around m_DM/20) and experience decoupling from the Universe plasma. WIMP candidates such as the neutralino <cit.>  and Kaluza-Klein particle <cit.> are found in theories such as the minimal supersymmetric Standard Model (MSSM) and universal extra dimensions (UED), respectively, and also in other extensions of the Standard Model (SM) such as singlet scalar <cit.> (or fermionic <cit.>) DM. In spite of their popularity,  WIMPs have not yet been  detected in  direct experiments. 

The other viable and well-motivated hypothesis to explain the DM problem is that there is such a feeble interaction that DM particles can never be abundant enough to thermalize. In this so-called freeze-in mechanism  <cit.>, feebly interacting massive particles (FIMPs)  have been  slowly produced in the early Universe through the collisions or decays of the bath particles. FIMP candidates are motivated in various extensions of the SM <cit.> and a well-known example which arises from  neutrino physics is the sterile neutrino <cit.>.  It is difficult to detect FIMP particles because of their small couplings with the SM. As for the indirect  searches, depending on the type of DM candidate, i.e. scalar <cit.>, fermion <cit.>, etc, some experiments have parameter space where they could survive but these are very borderline. For a study of the non-thermal properties of dark matter see Ref. <cit.>.

Although a lot of attention has been dedicated to single-particle DM models, some studies have considered DM models with the contribution of more particles in the observed DM density (multi-component DM  <cit.>). The simplest and the most common case is the union of the singlet scalar (fermionic) and the singlet fermionic (scalar) models which are employed in both freeze-in <cit.> and freeze-out <cit.> solutions (or intermediate cases <cit.>). Nevertheless, it remains a mystery whether DM is a single particle or multi-component.

In this paper, we analyze whether the freeze-in approach can
be properly used to produce the observed DM density in our Universe.  We choose a minimal two-component DM model in such a way that both of the components are FIMP particles. Following our hypothesis, we consider a singlet scalar and a Dirac fermion  where   an accidental symmetry guarantees their stabilities and a Higgs portal enables them to interact with the SM particles. The most striking feature of  our model is its simplicity,  as  the two candidates of DM particles do not  couple with each other and the model has separate overlaps with both the singlet scalar model <cit.> and the singlet fermionic model <cit.>.  In our work, all contributing processes to the relic density are assumed and supplementary phenomenological aspects  are also included.

Some promising possible signatures of FIMPs which are found to be most reliable in previous works are the γ-ray excess  observed from the Galactic center (GC) <cit.>, the X-ray line at 3.55 keV <cit.>, and DM self-interaction <cit.>. To generate the gamma ray excess, the fermionic component should have a pseudoscalar coupling to the mediator in the freeze-out regime <cit.>. The scalar component which does also couple to the SM Higgs directly <cit.>, should not feature large valued couplings. An X-ray signal with E_γ=3.55 keV from the XMM-Newton telescope and a similar signal at 3.52 keV from the Andromeda galaxy (M31) and Perseus Cluster could all  be interpreted by the decay <cit.> or the annihilation <cit.> of DM.  However, this requires a definite decay rate and annihilation cross section which is out of reach for our scenario. Therefore, we continue our probe relying only on the DM self-interaction. This non-gravitational interaction is a well-motivated indirect search  as it solves the tensions between observations and simulations of the small-scale structure of DM.    

Following the aforementioned setup, our paper is organized as follows. After introducing the construction of our model and identifying its parameter space in Section <ref>, we solve two independent Boltzmann equations in the following section (Section <ref>), in order to reach the observed relic density measured by the WMAP and Planck experiments <cit.>. In Section <ref>, we study the phenomenological implications for both direct and indirect experiments, and summarize our results in Section <ref>.                                             



§ TWO-COMPONENT FIMP DM



Beyond the SM, we employ three new fields to furnish our model: two scalars (χ and S) and one Dirac fermion (ψ), which are all  assumed to be singlet under the SM gauge groups. A discrete Z_2 symmetry is applied such that it reads the SM fields and  the S-scalar even, and the other two fields  (χ and ψ) odd. This symmetry guarantees the stability of both odd particles in a way that we do not have any terms involving both fields ψ and χ. In this way, the decays of odd particles to  one another  are prevented. Therefore, we can have two DM candidates in our setup by an accidental symmetry. 

The framework of  our model is constructed by:

    L    ⊃12(∂_μS)^2+12 (∂_μχ)^2+iψ̅∂ψ
       -m_ψψ̅ψ-g_sSψ̅ψ-g_pSψ̅γ^5ψ
       -V(H,S,χ),

where we introduced  the scalar and pseudoscalar interactions with the couplings g_s and g_p, respectively, and inserted the scalar interactions in the term  V(H,S,χ) as     

    V(H,S,χ)    =-μ_H^2H^†H+λ_H(H^†H)^2
       +μ_1S+12μ_S^2 S^2+13!α_s S^3+14!λ_S S^4+12m^2_0χχ^2+14!λ_χχ^4
       +λ_1S H^†H+ λ_2S^2 H^†H +λ_χ Hχ^2H^†H+λ_3Sχ^2+λ_4S^2χ^2.

After spontaneous symmetry breaking, the SU(2) Higgs doublet is parametrized as   

    H = 1/√(2)( [     0; v_H+h ]) ,
  
where v_H=246 GeV, but for the mediator we assume that it does not acquire a vacuum expectation value, i.e. <S>=0,  which minimalizes our model too. Now, due to the interaction terms in Eq. (<ref>), h and S mix with each other and form a mass matrix  with the following eigenstates

    h_1    = Ssinθ + hcosθ,
    
    h_2    = Scosθ - hsinθ,

and  the eigenvalues as

    m^2_h_1,h_2 = m^2_h+m^2_S/2±m^2_h-m^2_S/2√(1+y^2),     with   y= 2m^2_h,S/m^2_h-m^2_S,

where θ is the mixing angle between h_1 and h_2 such that

    tanθ = y/1+√(1+y^2).

According to the definition of the mixing angle θ, h_1 can be considered as the SM-like Higgs observed at the LHC with a mass of about 125 GeV.  In Eq. (<ref>), m_h=√(2λ_H)v_H, m_S=(λ_2v_H^2+μ_S^2)^1/2 and m_h,S=√(λ_1v_H).

Concerning our parameters, vacuum stability implies that the scalar potential in Eq. (<ref>) must be bounded from below. On the other hand, perturbativity does not allow the model parameters to be too large. Eventually, these theoretical conditions can be satisfied  if one has

    -2π/3   <λ_S, λ_χ<2π/3,
    
    -4π   <λ_2, λ_χ H<4π,
    
    -8π   <λ_4, g_s, g_p<8π,
    λ_2   +√(λ_Hλ_S)>0,
    λ_χ H   +√(λ_Hλ_χ)>0,
    
    2λ_4   +√(λ_Sλ_χ)>0,

and

    (√(2(λ_2+√(λ_Hλ_S))(λ_χ H+√(λ_Hλ_χ))(2λ_4+√(λ_Sλ_χ)))
    
    +   √(λ_Hλ_Sλ_χ)+λ_2√(λ_χ)+λ_χ H√(λ_S)+2λ_4√(λ_H))>0,

where λ_H is the quartic coupling of H. Extending the SM with the new fields ψ, χ and S embeds 19 parameters in addition to the SM ones.  They are  m_S, m_χ, m_ψ, g_s, g_p, μ_1, α_S, λ_S, λ_χ, λ_1, λ_2, λ_χ H, λ_3, λ_4, m_h_1, m_h_2, sinθ, v_H, and m_h.  However, due to the 8 model constraints, 11 independent parameters,  

    g_s, g_p, λ_3, λ_4, λ_χ H, m_ψ, m_χ, m_h_2,  sinθ,

remain for the relic abundance and for indirect searches, λ_H and λ_χ are required. Here, we take a moment to describe the eight constraints which appear in our work. Note that, after spontaneous symmetry breaking the scalar potential given in Eq. (<ref>) reads as V(h,S,χ). 
Therefore, it can be deduced from the potential that:

    1)        ∂ V∂ h|_h=S=χ=0=0 ⇒μ_H^2=λ_Hv_H^2,
    
    2)        ∂ V∂ S|_h=S=χ=0=0 ⇒μ_1=-λ_1v_H^22,
    
    3)         m_h^2=-μ_H^2+3λ_Hv_H^2=2λ_Hv_H^2,
    
    4)          m_S^2=μ_S^2+λ_2v_H^2,
    
    5)         m_χ^2=m^2_0χ+λ_χ Hv_H^2.
 
Also, the mixing between S and h produces the scalars h_1 and h_2 so one can conclude that 

    12m_h^2h^2+12m_S^2S^2+λ_1v_HSh=12m_h_1^2h_1^2+12m_h_2^2h_2^2.

Substituting h_1=Ssinθ + hcosθ and h_2=Scosθ - hsinθ (<ref>) and using the constraints (3)-(5), one obtains:

    6)      λ_H=m_h_1^2cos^2θ+m_h_2^2sin^2θ2v_H^2,
    
    7)      λ_2=m_h_1^2sin^2θ+m_h_2^2cos^2θ-μ_S^2v_H^2,
    
    8)      λ_1=m_h_1^2-m_h_2^22v_Hsin2θ.

These 8 constraints reduce the 19 free parameters in the model to the 11 independent parameters.
Also, the couplings α_s and λ_2 can be taken as zero without any ambiguities. However, we consider α_s≠ 0 for  future applications.  In the following, we will probe our model parameter space with experimental constraints coming from the relic density, direct and indirect detections.



§ DM DENSITY





 §.§ Boltzmann Equation



Since our model contains two DM candidates, its relic density has contributions from both fields ψ and χ. Therefore, we have to solve two Boltzmann equations for particles which will not reach equilibrium in the freeze-in mechanism where we follow the solution  in Ref. <cit.> (following Ref. <cit.>). The  time evolution of number density, dn_DM/dt, for the fermionic DM  is given by  

    dn_ψdt+3Hn_ψ   = Tπ^2∑_i=1^2m_h_i^2K_1(m_h_iT) Γ_h_i→ψ̅ψ
       +T32π^4∑_j=f,Z,W,h_1,h_2∫_4m_j^2^∞dsσ_jj→ψ̅ψ(s)(s-4m_j^2)√(s)K_1(√(s)T),
 
and for the scalar DM, it reads

    dn_χdt+3Hn_χ   = Tπ^2∑_i=1^2m_h_i^2K_1(m_h_iT)Γ_h_i→χχ
       +T32π^4∑_j=f,Z,W,h_1,h_2∫_4m_j^2^∞dsσ_jj→χχ(s)(s-4m_j^2)√(s)K_1(√(s)T).

Here H is the Hubble constant, K_1 is the modified Bessel function of order 1 and s is the  center of mass energy squared. All contributions to the DM relic density  are considered in the corresponding cross sections and decay widths in the two above equations. Our analytical results for the cross sections and decay widths are presented in the Appendix. The number density of DM particles is calculated as n_i=g_i(2π)^3∫ d^3pf_i <cit.> (with i=ψ, χ), where f_i is the phase space density of particle i with the g_i-internal spin degrees of freedom. As it is well-known from the freeze-in mechanism of production, the two DM candidates in the present model have negligible initial abundance (individually), thus we may set f_i = 0. Consequently, it can be derived from Eqs. (<ref>) and (<ref>)  that the process of DM conversion, i.e. χχ↔ψ̅ψ, does not contribute to the total relic abundance and is suppressed in our  next calculations. On the other hand, each of the DM candidates, independent of the other, can be produced or annihilated in the Universe.

By solving Eqs. (<ref>) and (<ref>), one can obtain the number densities (n_ψ, n_χ) scaled to the entropy of Universe ŝ, i.e. Y_ψ=n_ψ/ŝ and Y_χ=n_χ/ŝ, as    

    Y_ψ   = 14π^445M_pl1.66g^*_s(T)√(g^*_ρ)[2∑_i=1^2Γ_h_i→ψ̅ψm_h_i^2∫_T_Now^∞dTK_1(m_h_iT)T^5
       +∑_j=f,Z,W,h_1,h_2116π^2∫_T_Now^∞dT1T^5∫_4m_j^2^∞dsσ_jj→ψ̅ψ(s)(s-4m_j^2)√(s)K_1(√(s)T)],

and

    Y_χ   = 14π^445M_pl1.66g^*_s(T)√(g^*_ρ)[2∑_i=1^2Γ_h_i→χχm_h_i^2∫_T_Now^∞dTK_1(m_h_iT)T^5
       +∑_j=f,Z,W,h_1,h_2116π^2∫_T_Now^∞dT1T^5∫_4m_j^2^∞dsσ_jj→χχ(s)(s-4m_j^2)√(s)K_1(√(s)T)],
 
where M_pl is the Planck mass, and g^*_s and g^*_ρ are the effective numbers of degrees of freedom. 



 §.§ Relic Abundance



The most important constraint which should be satisfied in models describing  DM is the observed relic  density. As the Planck experiments have measured the current amount of DM <cit.>, our first experimental constraint is described as

    Ω_DMh^2=Ω_ψh^2+Ω_χh^2=0.1199±0.0027,
   
where h is the Hubble parameter scaled in units of 100 km/s.Mpc. Using the  yield calculated in the previous section (Eqs. (<ref>) and (<ref>)), we can obtain the relic density as  

    Ω_ih^2=2.742×10^-8(m_iGeV)Y_i(T_0),     i=ψ, χ.
 
First, we start with the scalar component of the model. The dependence of DM density is evaluated over the relevant parameters.  The predicted relic density of our model is best behaved at m_h_2=100 GeV and sinθ=0.01, with the required value of 5×10^-10 GeV for  mediator-scalar DM coupling λ_3. Two other couplings, λ_χ H and λ_4, are found to have major roles in controlling the relic density. By varying the singlet scalar DM mass (inspired by Eq. (<ref>)), we probe our parameter space in two classes: different values of λ_χ H and of λ_4 (see Fig. <ref>).       



The behavior of the relic density Ω_χ h^2 regarding different Higgs-scalar DM couplings, λ_χ H, is depicted  in the logarithmic scale in Fig. <ref>. In addition to the fixed values of relevant parameters λ_3, m_h_2 and sinθ, we have adopted λ_χ H as 10^-12, 5×10^-13 and 10^-13  as we  fixed λ_4=10^-13.   Through Fig. <ref>, it is obvious that the resonance occurs at m_χ∼ m_h_2/2. For masses below the resonance, the relic density of the scalar component increases linearly  in the log-scale, but for larger values (m_χ> m_h_2/2) it  seems that the relic density is independent of the mass m_χ. The difference between these two regions is  due to the process h_2→χχ which is allowed in the region below the resonance.         

The complementary analysis of the scalar component is plotted in Fig. <ref>, where we have chosen λ_4=10^-12, 5×10^-13 and 10^-13. Regarding the resonance  at m_χ∼ m_h_2/2, as in Fig. <ref>, it can be seen that for the region below the resonance, the relic density grows when the scalar mass m_χ increases. This part of the graph seems to be independent of the λ_4-value and it is enhanced by the h_2→χχ process. After a significant drop at m_χ∼ m_h_2/2, the relic density seems to be  independent of DM mass for the region above the resonance. It is mainly influenced by changing the quartic coupling λ_4.  Finishing our investigation of the scalar component, it should be noted that this analysis has a good overlap with a singlet scalar model <cit.>.        
We continue our investigation in parallel by turning our attention to the fermionic DM  in the logarithmic scale, too. As before, we consider two classes of variations defined by the effect of scalar (parameterized by g_s) and pseudoscalar (parameterized by g_p) interactions of ψ (see Fig. <ref>). We first look at the coupling g_s so  its best effects are formed for the values of 10^-8, 10^-9 and 10^-10 (Fig. <ref>). Similar to the scalar case, the resonance position occurs at m_ψ∼ m_h_2/2, so below this value one can observe the linear behavior of relic density which arises through the process h_2→ψψ. For m_ψ>m_h_2/2, the relic density changes by several order of magnitudes in a small interval of mass range. In this region, Ω_ψh^2 decreases when m_ψ increases. Note that the relic density is approximately independent of DM mass for  large  values of g_s. Also,  for small enough values of g_s, decreasing g_s does not significantly change the relic density.



We keep on probing our model parameter space by choosing  the appropriate values of g_p as 10^-9, 10^-10 and 5×10^-11.  In Fig. <ref>, we show the behavior of the relic density of the fermionic DM in terms of its mass. Here, there is a distinct point  which should be expressed. As is seen, for very small values of g_p, the relic density is approximately independent of mass for massive DM.  Similar to the scalar component, we can compare the fermionic component with models describing singlet fermionic DM like in Ref. <cit.>.
   


§ PHENOMENOLOGICAL IMPLICATIONS





 §.§ Direct Searches



In this section, we search for signals inspired from XENON100 <cit.> and LUX <cit.> in spin-independent elastic scattering of DM off nuclei. Our intended process includes the fundamental interaction of DM-quark which occurs via the t-channel mediated by scalars h_1 and h_2. Taking into account the contribution of each DM component and using fractions ξ_ψ=Ω_ψ/Ω_DM and ξ_χ=Ω_χ/Ω_DM, we investigate whether the model parameter space could be affected by the experimental results in this way.  To this end, we calculate the following cross sections,  


    σ_SI^ψ = ξ_ψ g_s^2μ_m^2sin^2θcos^2θ/π(1m_h_1^2-1m_h_2^2)^2λ_N^2,

and

    σ_SI^χ   = ξ_χμ_m^2/4π m_h_1^4m_h_2^4m_χ^2[2 v_H λ_χ H(m_h_1^2sin^2θ+m_h_2^2cos^2θ),
       +λ_3sin2θ(m_h_1^2-m_h_2^2)]^2λ_N^2,

where

    λ_N = m_N/v_H[∑_q=u,d,sf_q+29(1-∑_q=u,d,sf_q)]≈1.4×10^-3.

Here, the parameter μ_m=m_Nm_DMm_N+m_DM is the reduced mass of DM-nucleon and m_χ=(m_0χ^2+λ_χ Hv_H^2)^1/2 is the physical mass of the scalar DM. A cancellation effect <cit.> could occur when the two terms in Eq. (<ref>) cancel each other out, giving a suppressed cross section which is not appropriate for our consideration. Generally, as was  mentioned earlier, a necessary condition for our DM candidates to be nonthermal is  that they have extremely small couplings (g_s, λ_χ H) which would yield cross sections out of the sensitivity of the aforementioned experiments by their established values of order ∼10^-8-10^-12. Searching for other viable experiments, we consider the scattering of DM off free electrons in materials such as superconductors, semiconductors and graphene.  From Ref. <cit.> it is seen that, although, these  electron detectors are useful for light DM particles (O(MeV)),  the mediator mass should also be  of order O(MeV), which is in conflict with the current model including m_h_2=100 GeV.  For this reason, we are not able to probe such FIMP models directly. This outcome is consistent with the lack of direct experimental signals to date.   



 §.§ Indirect Searches
 


  §.§.§ Invisible Higgs Decays







Since the ATLAS and CMS have recorded the signature of the SM Higgs <cit.>, new searches have been prepared for DM phenomenology. This is done by considering the branching ratios of the Higgs, especially for decaying into light DM candidates,

    Br(h_1→Invisible)=Γ^inv(h_1→χχ)+Γ^inv(h_1→ψψ)Γ^SM+Γ^inv.
    
Regarding the experimental upper bound 0.23 for Br(h_1→Invisible) <cit.>, we see that the decays of Higgs to both DM ψ and χ are suppressed due to small couplings g_s and g_p in the former case, and small λ_4 and λ_χ H in the latter. However, another constraint comes from the decay of  our Higgs to h_2 (if kinematically possible, i.e. m_h_2 < m_h_1/2) whose decay rate could be calculated as 

    Γ(h_1→ h_2h_2)=e_3^28π m_h_1(1-4m_h_2^2m_h_1^2)^1/2Θ(m_h_1-2m_h_2),

where e_3 is the relevant vertex factor which is presented in the Appendix, see Eq. (<ref>).  From Eq. (<ref>), it can be seen that for sinθ=0.01 and m_h_1=125 GeV,
the result is sensitive to the choice of α_s.  Consequently, we investigate the behavior of the aforementioned decay rate regarding the mass of the Higgs h_1 and the relevant coupling α_s.  The parameter space of our Higgs sector is plotted in Figs. <ref> and <ref> where, respectively, we have calculated Br(h_1→ h_2h_2) as a function of m_h_2 and depicted the parameter space for the (α_s, m_h_2)-plane which is consistent with experimental measurements.

 


  §.§.§ DM Self-Interactions



Of the different DM models, the collision-less cold DM (CDM) paradigm has been successful in explaining the large scale structure of the Universe. However, there
are discrepancies between the CDM predictions and observations on smaller scales. The self-interacting DM (SIDM) paradigm has the potential to solve these issues (for a review of SIDM, see Ref. <cit.>). Although such interactions cannot be detected in experiments, we can infer bounds on σ_DM/m by evaluating the trajectory of DM in colliding
galaxy clusters <cit.>. An updated work <cit.> has considered a set of twelve galaxies and six clusters in order to cover different scales. Including the core sizes from dwarf to cluster (varying from 0.5 to 50 kpc), the aforementioned cross section is parametrized as

    σ_DM/m∼0.1-2 cm^2g^-1.

In this section, we analyze this constraint to see if it can put new limits on the parameter space of our model. 

The DM self-interaction in the present model includes the processes χχ→χχ, ψψ→ψψ, χχ→ψψ, ψψ→χχ and χψ→χψ. Except for χχ→χχ and ψψ→ψψ, the processes contain cross sections proportional to the coupling λ_4, which is very small in our work. Therefore, the specified processes do not contribute to this cosmological constraint. Concerning the processes χχ→χχ and ψψ→ψψ, we  start first with the scalar component which has been studied in a singlet FIMP scalar model in Ref. <cit.>.  Here, we just consider the contact interaction which is parameterized by the coupling λ_χ. Practically, we neglect the contributions from the s-channel mediated diagrams. This is due to the small couplings of the scalar DM with both the SM Higgs and the mediator and also due to the large masses which appear in the propagator. One way to vitalize the s-channel contribution might be through fine tuning by considering the scattering near resonance (similar to Ref. <cit.>). In this way, in the denominator of the propagator, m_h_2 should be tuned such that |m_χ-m_h_2|≪1 GeV. Considering the values of couplings needed  for the observed relic density,  this scenario  fails too. Therefore, following Ref. <cit.>, we obtain the  self-interaction cross section per mass m_χ, as

    σ_χm_χ=9λ_χ^22π m_χ^3,

where λ_χ is the quartic self-coupling of the scalar DM (see Eq. (<ref>)). Following the theoretical constraints in Eqs. (<ref>) and (<ref>), we obtain an experimental upper bound of about 0.1 GeV on the mass of the scalar DM, which is depicted in Fig. <ref>. Going back to Figs. <ref> and <ref>, we observe that this range of scalar mass can produce proper total relic density along with the contribution of the fermionic component.



Another significant point, which we would like to clarify in this work, is the self-interaction of a singlet fermionic FIMP DM. In general, we have two concerns. First, it should be noted that significant self-scattering at dwarf scales requires the mediator masses to be smaller than 100 MeV <cit.>. In fact, following Ref. <cit.>, the fermionic DM should satisfy the relation  (m_ψ /10  GeV) (m_h_2 / 100  MeV)^2∼1, with Yukawa interactions of strengths 10^-5 to 1, which is in contradiction with our fermionic DM coupling and mediator mass. 
Second, if the mediator couples to the SM through a Higgs portal, one should make sure that the mediator decays before the start of Big Bang nucleosynthesis (BBN), so the decay products do not affect the BBN. Eventually, we require a mediator with a lifetime ∼ 1 s. One way to alleviate the second constraint in DM models with extremely weak interactions is to open a new decay channel for the mediator so it can decay faster. This is done
in Refs. <cit.> and <cit.> by coupling the mediator to a light sterile neutrino (it should be noted that this new coupling does not affect the relic density). However, the first constraint (light mediator) is in conflict with our mediator of mass 100 GeV (and other usual two-component models). A promising solution seems to be to work at the resonance region  to minimize this mass constraint <cit.>. Due to the small coupling g_p (pseudoscalar interaction type), and the fact that there is an energy (velocity) dependent correction to the width in the resonance region (as explained in Ref. <cit.>), we conclude that the resonance DM self-interaction scenario does not work in fermionic FIMP models.               



§ CONCLUSIONS AND OUTLOOK



We have constructed a minimal two-component model to analyze the implications of multi-component DM in the Universe. Using the freeze-in mechanism, we calculated the relic abundance predicted by our model and compared it to the observed relic abundance of DM. We started our investigation by proposing two DM particles: a real scalar and a Dirac
fermion. Furthermore, a scalar mediator between the dark sector and the SM sector was added. The couplings for this interaction are assumed to be small as we are utilizing the freeze-in
mechanism. We solved two independent Boltzmann equations in order to obtain the observed relic density with  the contributions of both DM components.
It should be noted here that at the time of finishing this work, a new version of micrOMEGAs <cit.> was presented  which can compute the relic abundance of FIMP candidates. In the following, using theoretical constraints, we probed the model parameter space and compared our results with the relevant singlet models. Although it is difficult to probe FIMP particles, we looked for astrophysical probes, first considering direct detection.  We considered the scattering of DM particles off nuclei and free electrons. As we explained, it is impossible to see this direct signature for our FIMP model. 

In order to constrain the parameter space of our model, we also checked the limits from the invisible decay width of the Higgs.  Finally, we probed the self-interaction of DM
in this model. We used the bounds on non-gravitational interactions of DM in giant cluster collisions and constrained the mass of DM candidates in our model. In addition to the mentioned probes of DM, we can refer to the Big Bang nucleosynthesis  and cosmic microwave background  constraints regarding dark photons and dark Higgs <cit.>. These neutral bosons  mix with the SM photon (kinetically) and the SM Higgs, respectively, by accepting significant bounds on their coupling parameters. 

Two-component DM is a starting point for considering multi-component configurations where DM consists of various types of fundamental particles (scalar, fermion, vector and etc). The freeze-in framework is also a well-motivated approach which may be probed more extensively by future experiments.                           





§ ACKNOWLEDGMENT


We are particularly grateful to Yonit Hochberg for giving us insights into the direct probes, and we would like to thank Takashi Toma, Ian Shoemaker, Madhurima Pandey and Anirban Biswas for useful discussions.  







§ APPENDIX:  DM PRODUCTION CROSS SECTIONS AND DECAY RATES


Here, we present our calculation of the fermionic DM  production cross-sections which contribute to the relic density of our  model:


    (σ v_rel)_ff→ψψ=ssin^22θ/32πN_c(m_f/v_H)^2(1-4m_f^2/s) F



    (σ v_rel)_ZZ→ψψ=sin^22θ/36π(m_Z^2/v_H)^2(2+(s-2m_Z^2)^2/4m_Z^4)F



    (σ v_rel)_W^+W^-→ψψ=sin^22θ/36π(m_W^2/v_H)^2(2+(s-2m_W^2)^2/4m_W^4)F

where


    F   =   (1-4m_ψ^2/s)^3/2(g_ψ^2+sg_p^2/s-4m_ψ^2){1/(s-m_h_1^2)^2+m_h_1^2Γ_h_1^2+
          1/(s-m_h_2^2)^2+m_h_2^2Γ_h_2^2-2(s-m_h_1^2)(s-m_h_2^2)+2m_h_1m_h_2Γ_h_1Γ_h_2/[(s-m_h_1^2)^2+m_h_1^2Γ_h_1^2][(s-m_h_2^2)^2+m_h_2^2Γ_h_2^2]}.



    (σ v_rel)_ h_1h_1→ψψ   =   1/4π s(1-4m_ψ^2/s)^1/2{(g_p^2s+g_ψ^2(s-4m_ψ^2))
          ×[e_1^2sin^2θ/(s-m_h_1^2)^2+m_h_1^2Γ_h_1^2+e_2^2cos^2θ/(s-m_h_2^2)^2+m_h_2^2Γ_h_2^2
          +e_1e_2sin2θ((s-m_h_1^2)(s-m_h_2^2)+m_h_1m_h_2Γ_h_1Γ_h_2)/((s-m_h_1^2)^2+m_h_1^2Γ_h_1^2)((s-m_h_2^2)^2+m_h_2^2Γ_h_2^2)]
          -8m_ψ g_ψsin^2θ[e_1sinθ(s-m_h_1^2)/(s-m_h_1^2)^2+m_h_1^2Γ_h_1^2+e_2cosθ(s-m_h_2^2)/(s-m_h_2^2)^2+m_h_2^2Γ_h_2^2]
          ×[arctanh(x_1)[g_p^2(2m_h_1^2-3s)-g_ψ^2(2m_h_1^2-8m_ψ^2+s)]/x_1(s-2m_h_1)+g_p^2-g_ψ^2]
          -4sin^4θ[g_p^4+g_ψ^4-2g_p^2g_ψ^2+(m_h_1^2(g_p^2-g_ψ^2)+4g_ψ^2m_ψ^2)^2-8g_p^2g_ψ^2m_ψ^2s/m_h_1^4+m_ψ^2(s-4m_h_1^2)
          -arctanh(x_1)/x_1(s-2m_h_1)^2(g_p^4[6m_h_1^4+s(s-4m_h_1^2)]
          -2g_p^2g_ψ^2[2m_h_1^2(3m_h_1^2-4m_ψ^2)+s(s-4m_h_1^2)]
          +g_ψ^4[2m_h_1^2(3m_h_1^2-8m_ψ^2)+s(s-4m_h_1^2)+16m_ψ^2(s-2m_ψ^2)])]},



    (σ v_rel)_ h_2h_2→ψψ   =   1/4π s(1-4m_ψ^2/s)^1/2×{(g_p^2s+g_ψ^2(s-4m_ψ^2))
          ×[e_3^2sin^2θ/(s-m_h_1^2)^2+m_h_1^2Γ_h_1^2+e_4^2cos^2θ/(s-m_h_2^2)^2+m_h_2^2Γ_h_2^2
          +e_3e_4sin2θ[(s-m_h_1^2)(s-m_h_2^2)+m_h_1m_h_2Γ_h_1Γ_h_2]/[(s-m_h_1^2)^2+m_h_1^2Γ_h_1^2][(s-m_h_2^2)^2+m_h_2^2Γ_h_2^2]]
          -8m_ψ g_ψcos^2θ[e_3sinθ(s-m_h_1^2)/(s-m_h_1^2)^2+m_h_1^2Γ_h_1^2+e_4cosθ(s-m_h_2^2)/(s-m_h_2^2)^2+m_h_2^2Γ_h_2^2]
          ×[arctanh(x_2)[g_p^2(2m_h_2^2-3s)-g_ψ^2(2m_h_2^2-8m_ψ^2+s)]/x_2(s-2m_h_2)+g_p^2-g_ψ^2]
          -4cos^4θ[g_p^4+g_ψ^4-2g_p^2g_ψ^2+(m_h_2^2(g_p^2-g_ψ^2)+4g_ψ^2m_ψ^2)^2-8g_p^2g_ψ^2m_ψ^2s/2m_h_2^4+2m_ψ^2(s-4m_h_2^2)
          -arctanh(x_2)/x_2(s-2m_h_2)^2(g_p^4(6m_h_2^4+s(s-4m_h_2^2))
          -2g_p^2g_ψ^2(2m_h_2^2(3m_h_2^2-4m_ψ^2)+s(s-4m_h_2^2))
          +g_ψ^4(2m_h_2^2(3m_h_2^2-8m_ψ^2)+s(s-4m_h_2^2)+16m_ψ^2(s-2m_ψ^2)))]},



    (σ v_rel)_ h_1h_2→ψψ   =   1/4π s(1-4m_ψ^2/s)^1/2{[g_p^2s+g_ψ^2(s-4m_ψ^2)]
          ×[e_2^2sin^2θ/(s-m_h_1^2)^2+m_h_1^2Γ_h_1^2+e_3^2cos^2θ/(s-m_h_2^2)^2+m_h_2^2Γ_h_2^2
          +e_2e_3sin2θ((s-m_h_1^2)(s-m_h_2^2)+m_h_1m_h_2Γ_h_1Γ_h_2)/[(s-m_h_1^2)^2+m_h_1^2Γ_h_1^2][(s-m_h_2^2)^2+m_h_2^2Γ_h_2^2]]
          -4m_ψ g_ψsin2θ[e_2sinθ(s-m_h_1^2)/(s-m_h_1^2)^2+m_h_1^2Γ_h_1^2+e_3cosθ(s-m_h_2^2)/(s-m_h_2^2)^2+m_h_2^2Γ_h_2^2]
          ×[arctanh(x_i)/x_i(s-2m_h_i^2)[g_p^2(m_h_1^2+m_h_2^2-3s)-g_ψ^2(m_h_1^2+m_h_2^2-8m_ψ^2+s)]+g_p^2-g_ψ^2]
          -sin^22θ[1/2m_h_i^4+2m_ψ^2(s-4m_h_i^2)(-g_p^4[2m_h_i^4+m_h_i^2(m_h_j^2-8m_ψ^2)+2m_ψ^2s]
          +2g_p^2g_ψ^2(2m_h_i^4+m_h_i^2(m_h_j^2-10m_ψ^2)-2m_ψ^2(m_h_j^2-3s)))
          -g_ψ^4[(m_h_i^2-4m_ψ^2)(2m_h_i^2+m_h_j^2-4m_ψ^2)-2m_ψ^2s]
          +arctanh(x_i)/x_i(s-2m_h_i)^2(g_p^4[3m_h_i^2+m_h_i^2
    (4m_h_j^2-3s)-m_h_j^4-m_h_i^2s+s^2]
          +2g_p^2g_ψ^2[-3m_h_i^2+m_h_i^2(
    3(4m_ψ^2+s)-4m_h_j^4)+m_h_j^4+m_h_j^2(s-4m_ψ^2)-s^2]
          +g_ψ^4[3m_h_i^2+m_h_i^2(
    -3(8m_ψ^2+s)+4m_h_j^2)-m_h_j^4-m_h_j^2(s-8m_ψ^2)+s^2
          +16m_ψ^2(s-2m_ψ^2)])]},


where the auxiliary parameters and coupling constants have the following expressions   


    x_i=(s-4m_ψ^2)^1/2 (s-4m_h_i^2)^1/2/(s-2m_h_i^2),      with   i, j=1,2  and i≠ j,



    e_1=-sin^2θ( α_ssinθ+6v_Hλ_2cosθ)
    -3cos^2θ(λ_1sinθ+2v_Hλ_Hcosθ),



    e_2   =   sin2θ(-1/2α_ssinθ+3λ_Hv_Hcosθ)
          -λ_1cosθ(1-3sin^2θ)+2v_Hλ_2sinθ(1-3cos^2θ),



    e_3   =   sin2θ(-1/2α_scosθ-3λ_Hv_Hsinθ)
          -λ_1sinθ(1-3cos^2θ)-2v_Hλ_2cosθ(1-3sin^2θ),



    e_4=-cos^2θ( α_scosθ-6v_Hλ_2sinθ)
    -3sin^2θ(λ_1cosθ-2v_Hλ_Hsinθ).


The scalar component will account for the DM phenomenology by the following annihilation cross sections:


    (σ v_rel)_f̅f→χχ   =   1/16π(m_f/v_H)^2(1-4m_f^2/s)^3/2N_cR



    (σ v_rel)_ZZ→χχ   =   1/36π s(m_Z^2/v_H)^2(2+(s-2m_Z^2)^2/4m_Z^4)(1-4m_Z^2/s)^1/2R,



    (σ v_rel)_W^+W^-→χχ   =   1/36π s(m_W^2/v_H)^2(2+(s-2m_W^2)^2/4m_W^4)(1-4m_W^2/s)^1/2R


where


    R   =   e_5^2cos^2θ/(s-m_h_1^2)^2+m_h_1^2Γ_h_1^2+e_6^2sin^2θ/(s-m_h_2^2)^2+m_h_2^2Γ_h_2^2
          -e_5e_6sin2θ[(s-m_h_1^2)(s-m_h_2^2)+m_h_1m_h_2Γ_h_1Γ_h_2]/[(s-m_h_1^2)^2+m_h_1^2Γ_h_1^2][(s-m_h_2^2)^2+m_h_2^2Γ_h_2^2].



    (σ v_rel)_h_1h_1 →χχ   =   1/16π s(1-4m_h_1^2/s)^1/2
          ×{[e_1^2e_5^2/(s-m_h_1^2)^2+m_h_1^2Γ_h_1^2+e_2^2e_6^2/(s-m_h_2^2)^2+m_h_2^2Γ_h_2^2
          +2e_1e_2e_5e_6[(s-m_h_1^2)(s-m_h_2^2)+m_h_1m_h_2Γ_h_1Γ_h_2]/[(s-m_h_1^2)^2+m_h_1^2Γ_h_1^2][(s-m_h_2^2)^2+m_h_2^2Γ_h_2^2]]
          -(8e_5^2s-2m_h_1^2F(y_1)+2e_7)[e_1e_5(s-m_h_1^2)/(s-m_h_1^2)^2+m_h_1^2Γ_h_1^2
          +e_2e_6(s-m_h_2^2)/(s-m_h_2^2)^2+m_h_2^2Γ_h_2^2]
          +8e_5^2s-2m_h_1^2F(y_1)[e_5^2s-2m_h_1^2[1F(y_1)(1-y_1^2)+1]+e_7]
    +e_7^2},



    (σ v_rel)_h_2h_2 →χχ   =   1/16π s(1-4m_h_2^2/s)^1/2
          ×{[e_3^2e_5^2/(s-m_h_1^2)^2+m_h_1^2Γ_h_1^2+e_4^2e_6^2/(s-m_h_2^2)^2+m_h_2^2Γ_h_2^2
          +2e_3e_4e_5e_6[(s-m_h_1^2)(s-m_h_2^2)+m_h_1m_h_2Γ_h_1Γ_h_2]/[(s-m_h_1^2)^2+m_h_1^2Γ_h_1^2][(s-m_h_2^2)^2+m_h_2^2Γ_h_2^2]]
          -(8e_6^2s-2m_h_2^2F(y_2)+2e_8)[e_3e_5(s-m_h_1^2)/(s-m_h_1^2)^2+m_h_1^2Γ_h_1^2
          +e_4e_6(s-m_h_2^2)/(s-m_h_2^2)^2+m_h_2^2Γ_h_2^2]
          +8e_6^2s-2m_h_2^2F(y_2)[e_6^2s-2m_h_2^2(1F(y_2)(1-y_2^2)+1)+e_8]
    +e_8^2},



    (σ v_rel)_ h_1h_2→χχ   =   1/16π s(1-4m_h_i^2/s)^1/2
          ×{[e_2^2e_5^2/(s-m_h_1^2)^2+m_h_1^2Γ_h_1^2+e_3^2e_6^2/(s-m_h_2^2)^2+m_h_2^2Γ_h_2^2
          +2e_2e_3e_5e_6[(s-m_h_1^2)(s-m_h_2^2)+m_h_1m_h_2Γ_h_1Γ_h_2]/[(s-m_h_1^2)^2+m_h_1^2Γ_h_1^2][(s-m_h_2^2)^2+m_h_2^2Γ_h_2^2]]
          -(8e_5e_6s-2m_h_i^2F(y_i)+2e_9)[e_2e_5(s-m_h_1^2)/(s-m_h_1^2)^2+m_h_1^2Γ_h_1^2
          +e_3e_6(s-m_h_2^2)/(s-m_h_2^2)^2+m_h_2^2Γ_h_2^2]
          +8e_5e_6s-2m_h_i^2F(y_i)[e_5e_6s-2m_h_i^2(1F(y_i)(1-y_i^2)+1)+e_9]
    +e_9^2},


where we have employed parameter y_i and function F(y_i) as


    y_i=(s-4m_χ^2)^1/2 (s-4m_h_i^2)^1/2/(s-2m_h_i^2)
          F(y_i)=1/y_iarctanh(y_i)       with   i=1,2 ,


and also for coupling constants, we have the following parameters:


    e_5=-[v_Hλ_χ Hcosθ+λ_3sinθ]×2!,



    e_6=[v_Hλ_χ Hsinθ-λ_3cosθ]×2!,



    e_7=-[1/2λ_χ Hcos^2θ+λ_4sin^2θ]×4!,



    e_8=-[1/2λ_χ Hsin^2θ+λ_4cos^2θ]×4!,



    e_9=[(1/2λ_χ H-λ_4)sin2θ]×2!.


Finally, the decay rates of scalars h_i (with i=1,2) into fermionic and scalar DM particles are given as: 


    Γ(h_i→ψ̅ψ)=S_i^2θ/8π(1-4m_ψ^2/m_h_i^2)^3/2[g_ψ^2+m_h_i^2g_p^2/m_h_i^2-4m_ψ^2],



    Γ(h_1→χχ)=e_5^2/32π m_h_1(1-4m_χ^2/m_h_1^2)^1/2,



    Γ(h_2→χχ)=e_6^2/32π m_h_2(1-4m_χ^2/m_h_2^2)^1/2,

where we have defined S_1θ=sinθ and S_2θ=cosθ.




 
10


Gondolo:1990dk
P. Gondolo, and G. Gelmini,

Nucl. Phys. B, 360:145 (1991).




Srednickii
M. Srednicki, R. Watkins, and K. A. Olive,

Nucl. Phys. B, 310:693 (1988).




Chiu
H. Y. Chiu,

Phys. Rev. Lett,  17:712 (1966).





Jungman:1995df
G. Jungman, M. Kamionkowski, and K. Griest,

Phys. Rept,  267:195 (1996).






Cheng:2002ej
H. C. Cheng, J. L. Feng, and K. T. Matchev,

Phys. Rev. Lett,  89:211301 (2002).






Servant:2002aq
G. Servant, and T. M. P. Tait,

Nucl. Phys. B, 650:391 (2003).







Silveira:1985rk
V. Silveira, and A. Zee,

Phys. Lett. B,  161:136 (1985).




McDonald
J. McDonald,

Phys. Rev. D, 50: 3637 (1994).





Burgess
C. P. Burgess, M. Pospelov, and T. ter Veldhuis,

Nucl. Phys. B, 619:709 (2001).




Barger
V. Barger, P. Langacker, M. McCaskey et al,

Phys. Rev. D, 77:035005 (2008).






Kim:2008pp
Y. G. Kim, K. Y. Lee, and S. Shin,

JHEP, 0805:100 (2008).





Ettefaghi
M. M. Ettefaghi, and R. Moazzemi,

JCAP, 1302:048 (2013).







Fairbairn
M. Fairbairn, and R. Hogan,

JHEP, 1309:022 (2013).







McDonald:2001vt
J. McDonald,

Phys. Rev. Lett,  88:091304 (2002).








Hall:2009bx
L. J. Hall, K. Jedamzik, J. March-Russell, and S. M. West,

JHEP, 1003:080 (2010).






Bernal:2017kxu
N. Bernal, M. Heikinheimo, T. Tenkanen, K. Tuominen and V. Vaskonen,

Int. J. Mod. Phys. A 32 (2017) no.27,  1730023
doi:10.1142/S0217751X1730023X
[arXiv:1706.07442 [hep-ph]].





Yaguna:2011qn
 C. E. Yaguna,

JHEP, 1108:060 (2011).







Klasen:2013ypa
 M. Klasen, and C. E. Yaguna,

JCAP, 1311:039 (2013).







Ayazi:2015jij
S. Yaser Ayazi, S. M. Firouzabadi, and S. P. Zakeri,

J. Phys. G, 43 (9):095006 (2016).








Merle:2014xpa
A. Merle, and A. Schneider,

Phys. Lett. B, 749:283 (2015).




Merle
A. Merle, and M. Totzauer,

JCAP, 1506:011 (2015).



Shakya
B. Shakya,

Mod. Phys. Lett. A, 31 (06):1630005 (2016).




Kang:2014cia
Z. Kang,

Eur. Phys. J. C 75 (2015) no.10,  471
doi:10.1140/epjc/s10052-015-3702-4
[arXiv:1411.2773 [hep-ph]].



Biswas
A. Biswas, and A. Gupta,

JCAP, 1609:044 (2016).







Pandey:2017quk
M. Pandey, D. Majumdar, and K. P. Modak,

arXiv: hepph/1709.05955.



Dev:2013yza
P. S. Bhupal Dev, A. Mazumdar and S. Qutub,

Front. in Phys.  2 (2014) 26
doi:10.3389/fphy.2014.00026
[arXiv:1311.5297 [hep-ph]].




Profumo:2009tb
S. Profumo, K. Sigurdson, and L. Ubaldi,

JCAP, 0912:016 (2009).



Gelm1
G. B. Gelmini,

Nucl. Phys. Proc. Suppl,  138:32 (2005).



Gelm2
G. Duda, G. Gelmini, P. Gondolo et al,

Phys. Rev. D, 67:023505 (2003).



Dudaa
G. Duda, G. Gelmini, and P. Gondolo,

Phys. Lett. B, 529:187 (2002).





Herrero-Garcia:2017vrl
J. Herrero-Garcia, A. Scaffidi, M. White et al,

JCAP, 1711:021  (2017).





Biswas:2015sva
 A. Biswas, D. Majumdar, and P. Roy,

JHEP, 1504:065 (2015).





Esch:2014jpa
S. Esch, M. Klasen, and C. E. Yaguna,

JHEP, 1409:108 (2014).







Bhattacharya:2013hva
 S. Bhattacharya, A. Drozd, B. Grzadkowski et al,

JHEP, 1310:158 (2013).







Biswas:2013nn
 A. Biswas, D. Majumdar, A. Sil et al,

JCAP, 1312:049 (2013).






Chialva:2012rq
D. Chialva, P. S. B. Dev and A. Mazumdar,

Phys. Rev. D 87 (2013) no.6,  063522
doi:10.1103/PhysRevD.87.063522
[arXiv:1211.0250 [hep-ph]].




DuttaBanik:2016jzv
A. Dutta Banik, M. Pandey, D. Majumdar et al,

Eur. Phys. J. C, 77 (10):657 (2017).






Modak:2014vva
K. P. Modak,

JHEP, 1503:064 (2015).






Babu:2014pxa
K. S. Babu, and R. N. Mohapatra,

Phys. Rev. D, 89:115011 (2014).





Ade:2015dga
 P. A. R. Ade  et al (Planck Collaboration),

Astron. Astrophys,  571:A31 (2014).





Kolb:1990vq
E. W. Kolb and M. S. Turner,

Front. Phys.  69 (1990) 1.




Edsjo:1997bg 
J. Edsjo and P. Gondolo,

Phys. Rev. D 56, 1879 (1997).






Bernal:2017kxu
N. Bernal, M. Heikinheimo, T. Tenkanen, K. Tuominen and V. Vaskonen,

Int. J. Mod. Phys. A 32 (2017) no.27,  1730023.








Aprile:2012nq
E. Aprile et al (XENON100 Collaboration),

Phys. Rev. Lett,  109:181301 (2012).





Akerib:2013tjd
D. S. Akerib et al (LUX Collaboration),

Phys. Rev. Lett,  112:091303 (2014).





Hochberg:2015pha
Y. Hochberg, Y. Zhao, and K. M. Zurek,

Phys. Rev. Lett,  116 (1):011301 (2016).



Hoch1
Y. Hochberg, M. Pyle, Y. Zhao et al,

JHEP, 1608:057 (2016).



Hoch2
Y. Hochberg, Y. Kahn, M. Lisanti et al,

Phys. Lett. B, 772:239 (2017).





Aad:2012tfa
G. Aad  et al (ATLAS Collaboration),

Phys. Lett. B, 716:1 (2012).



Chatr
S. Chatrchyan et al (CMS Collaboration),

Phys. Lett. B, 716:30 (2012).





Belanger:2013kya
G. Belanger, B. Dumont, U. Ellwanger et al, 

Phys. Lett. B, 723:340 (2013).





Tulin:2017ara
S. Tulin, and H. B. Yu,

arXiv: hepph/1705.02358.




Clowe:2003tk
D. Clowe, A. Gonzalez, and M. Markevitch,

Astrophys. J,  604:596 (2004).





Randall:2007ph
S. W. Randall, M. Markevitch, D. Clowe et al,

Astrophys. J,  679:1173 (2008).





Kaplinghat:2015aga
M. Kaplinghat, S. Tulin and H. B. Yu,

Phys. Rev. Lett.  116 (2016) no.4,  041302.








Campbell:2015fra
R. Campbell, S. Godfrey, H. E. Logan et al, 

Phys. Rev. D, 92 (5):055031 (2015).





Tulin:2013teo
S. Tulin, H. B. Yu, and K. M. Zurek,

Phys. Rev. D, 87 (11):115007 (2013).





Kainulainen:2015sva
K. Kainulainen, K. Tuominen, and V. Vaskonen,

Phys. Rev. D, 93 (1):015016 (2016)
Erratum: [Phys. Rev. D, 95 (7):079901 (2017)].






Kouvaris:2014uoa
C. Kouvaris, I. M. Shoemaker, and K. Tuominen,

Phys. Rev. D, 91 (4):043519 (2015).






Duch:2017nbe
M. Duch, and B. Grzadkowski,

JHEP, 1709:159 (2017).





Belanger:2018ccd
G. Bélanger, F. Boudjema, A. Goudelis et al,

arXiv: hepph/1801.03509.




Fradette:2014sza
A. Fradette, M. Pospelov, J. Pradler and A. Ritz,

Phys. Rev. D 90 (2014) no.3,  035022.






Berger:2016vxi
J. Berger, K. Jedamzik and D. G. E. Walker,

JCAP 1611 (2016) 032.






