Light Fermionic WIMP Dark Matter with Light Scalar Mediator

A light fermionic weakly interacting massive particle (WIMP) dark matter is investigated by studying its minimal renormalizable model, where it requires a scalar mediator to have an interaction between the WIMP and standard model particles. We perform a comprehensive likelihood analysis of the model involving all robust constraints obtained so far and that will be obtained in the near future with paying particular attention to properly take the kinematically equilibrium condition into account. It is shown that near-future experiments and observations such as low-mass direct dark matter detections, flavor experiments and CMB observations play important roles to test the model; however, wide parameter region will still remain survived even if no WIMP and mediator signals are detected there. We also show that precise Higgs boson measurements at future lepton colliders will play a significant role to test this remaining region.


Introduction
In past decades, people tried to develop particle physics based on the electroweak naturalness [1,2], namely how the electroweak scale should be naturally explained. Many new physics scenarios such as supersymmetry, extra-dimension and composite Higgs have been proposed in this context, and those have been and still are being tested by various experiments including Large Hadron Collider (LHC). However, people recently start doubting this guiding principle, for new physics signals predicted by the scenarios have not been detected at all. In other words, though the electroweak scale should be naturally explained, it may be achieved by some other mechanisms (or ideas) which are totally different from what we have thought about so far. Particle physicists are seeking new mechanisms based on this consideration, but none of candidate models can successfully explain electroweak naturalness yet. Under this circumstance, one starts taking another strategy: developing particle physics by solving the dark matter problem. Once the nature of dark matter is clarified, it may launch out into the exploration of new physics beyond the standard model (SM).
The thermal dark matter, often called the weakly interacting massive particle (WIMP), is known to be one of influential dark matter candidates among others, for the dark matter abundance observed today is naturally explained by the so-called freeze-out mechanism [3,4], which can also successfully explain the history of the Big Band Nucleosynthesis (BBN) and recombination in the early universe. Though the WIMP is in general predicted to be in the mass range between (1) MeV [5,6] and (100) TeV [7][8][9][10][11][12], those with the mass around the electroweak scale have been intensively studied because of a possible connection to new physics models for the electroweak naturalness. However, present negative experimental results, not only from the LHC experiment but also from direct dark matter detection experiments, start eroding the parameter space of the WIMP with the electroweak mass. Thus, it motivates us to consider other WIMPs with a lighter ( (10) GeV) or heavier ( (1) TeV) mass. We focus on the former case in this paper.
Light WIMP must be singlet under the SM gauge group, otherwise it would be discovered already. Concerning the spin of the WIMP, we take one-half, namely a light fermionic WIMP in this paper. #1 In the minimal (renormalizable) model to describe such a light fermionic WIMP, a new additional particle called the mediator must be introduced to have an interaction between the WIMP and SM particles. In addition, such a mediator is required to be as light as the WIMP to explain the dark matter abundance observed today, to be singlet under the SM gauge group to avoid constraints from the current performed collider experiments, to be bosonic being consistent with the Lorentz symmetry, and to be even under the Z 2 symmetry in order to make the WIMP stable. Such a light fermionic dark matter with a light bosonic mediator recently receives many attentions, as it has a potential to have a large and velocity-dependent scattering cross section between WIMPs and solve the so-called small scale crisis of the universe [20][21][22]. Among two possibilities of the bosonic mediator, either a scalar or a vector [23][24][25][26][27], we take the scalar one in this paper. #2 To investigate the present status and future prospects of a light fermionic WIMP with a light scalar mediator, we perform a comprehensive analysis of its minimal (renormalizable) model [13,14,[30][31][32][33][34][35][36][37][38][39][40], where our likelihood involves all robust constraints obtained #1 Light scalar WIMP in its minimal model (Higgs-portal dark matter) is already excluded by the constraint from the invisible Higgs decay [13][14][15][16][17]. Light scalar WIMP still survives in next-to-minimal models [18,19]. #2 A careful model-building is required to have a light fermionic WIMP with a light vector mediator [28], because such a WIMP has a s-wave annihilation and tends to be excluded by the cosmic microwave background (CMB) observation. Such a constraint can be avoided for the scalar mediator case [29], as seen in section 3. so far and those will be obtained in the near future (if no WIMP and mediator signals are detected) from particle physics experiments as well as cosmological and astrophysical observations. We carefully involve a kinematical equilibrium condition assuming that the freezeout (chemical decoupling) of the WIMP occurs when it is in kinematically equilibrium with SM particles. We pay particular attention to a possible case that the light WIMP can be in the kinematical equilibrium via existent mediators at the freeze-out even if the WIMP does not have an interaction to SM particles with enough magnitude [41,42]. We find that a very wide parameter region is surviving at present in the certain mass region of the WIMP. We also show quantitatively how near-future experiments and observations such as low-mass direct dark matter detections, flavor experiments and CMB observations play important roles, by comparing the results of analyses for the present status and future prospects of the model. Moreover, we see that a wide parameter region will still remain even if neither WIMP nor mediator signal is detected in the near future, and show that precise Higgs boson measurements at future lepton colliders will play a significant role to test the region. Such comprehensive analysis by the global scanning is a recent trend in WIMP studies [43][44][45][46][47][48].
The rest of the paper is organized as follow. In section 2, we show our setup, the minimal renormalizable model to describe a light fermionic WIMP with a scalar mediator. We will give all interactions predicted by the model and discuss physics of the mediator. All constraints that we have involved in our likelihood are discussed in section 3. #3 Results of our likelihood analysis are given and discussed in section 4, including several implications of the result to (near) future projects for the WIMP search. Section 5 is devoted to the summary of our discussion. There are several appendices at the end of this paper, where preselection criteria we have involved in our analysis (Appendix A), the kinematical equilibrium condition (Appendix B) and results of our analysis (Appendix C) are shown in details.

The minimal model 2.1 Lagrangian
Two general properties of the WIMP field are considered before constructing a simplified model: the spin and weak-isospin(s) of the field(s). #4 Since we are interested in a light fermionic WIMP with its mass of (1) GeV or less, we focus on a Majorana WIMP, namely the simplest one among spin half WIMPs. On the other hand, the weak-isospin of the WIMP must be fixed to be zero, because a light WIMP carrying a non-zero weak-charge is excluded by collider experiments (LEP, etc.) performed so far. Such a singlet Majorana WIMP is wellmotivated by e.g. a neutralino (Bino, Singlino, etc.) in supersymmetric models.
An additional new mediator has to be introduced, otherwise the WIMP cannot have any renormalizable interaction with SM particles due to SM gauge symmetry, Lorentz symmetry and Z 2 symmetry making the WIMP stable. The mediator is required to be as light as the WIMP to satisfy the relic abundance condition, as seen in following sections. If we consider a Z 2 -odd mediator, it must be charged under the SM gauge interactions due to the SM gauge symmetry and such a light charged particle is also excluded by the collider experiments [50,51]. Hence, the mediator must be even under the Z 2 symmetry. Moreover, the mediator is #3 For readers who are not very much interested in the detail of the constraints and want to see the results of our analysis quickly, please skip this section (section 3) and go to the next section (section 4) directly. #4 Several dark matter fields with a different weak-charge are introduced for the well-tempered WIMP [49].
either a scalar or a vector boson, because it must have a renormalizable interaction with the WIMP (the WIMP-WIMP-mediator interaction). We consider the case of the scalar mediator in this paper, and leave the vector mediator case for future work.
We assume that the mediator is described by a real singlet from the viewpoint of minimality. Then, the Lagrangian involving all possible renormalizable interactions of the singlet Majorana WIMP χ, the mediator Φ and SM particles is given as follows [52]: with SM and H being the SM Lagrangian and the SM Higgs doublet, respectively. In order to make the WIMP stable as mentioned above, a Z 2 symmetry is imposed, where χ is odd under the symmetry but other particles are charged even. The scalar potential of the model is the potential of the SM doublet H involved in SM . Its explicit form is written as follows: Here λ i s are dimensionless coupling constants of quartic interactions, while others (µ i s and A ΦH ) are mass dimension one coupling constants for cubic and quadratic interactions.
We take v H = (−2µ 2 H /λ H ) 1/2 246 GeV and v Φ as vacuum expectation values of H and Φ. Taking the unitary gauge, the fields are expressed as H = [0, (v H +h )/ 2] T and Φ = v Φ +φ , where v Φ can be fixed to be zero without a loss of generality. Mass eigenstates of the scalars are then obtained by diagonalizing the quadratic terms of the potential, where The diagonalization matrix is described by the mixing angle θ , which controls the strength of interactions between φ and SM fermions or φ and SM gauge bosons. This angle is defined through the relation: Shown in following sections, the mixing angle sin θ is severely constrained, and hence its absolute value is phenomenologically required to be much smaller than one, namely |θ | 1. As a result, the mass eigenstates and the mixing angle are expressed as . (5) Note that the solution of the second equation for the angle θ has a multi-fold ambiguity within the domain of definition, −π ≤ θ ≤ π. However, we do not have to worry about this ambiguity because the angle is phenomenologically required to be suppressed as |θ | 1.

Interactions
Because of the mixing between the singlet Φ and the SM doublet H, the model predicts various interactions. First, interactions between the Higgs boson h and SM fermions & SM gauge bosons are suppressed slightly by cos θ compared to the SM prediction. On the other hand, the mixing introduces interactions between the mediator φ and the SM particles, where φ behaves like a light Higgs boson with the couplings suppressed by sin θ compared to the SM prediction. Both scalars have interactions with the WIMP as follows: Other interactions among SM fermions and SM gauge bosons in SM are not changed.
We next consider interactions among the scalars h and φ. From the scalar potential, four kinds of triple scalar interactions are obtained. Their explicit forms are given by θ + 6λ ΦH c 2 θ s 2 θ + λ Φ c 4 θ . As seen from all interactions discussed in this subsection, when the mixing angle θ is very suppressed, the scalar h becomes almost the SM Higgs boson, and it couples to the WIMP weakly. On the other hand, φ couples to the WIMP without any suppression, while interacts with SM particles very weakly. #5 As a result, the WIMP χ and the mediator φ form a dark sector coupling to the SM sector weakly through the mixing between H and Φ. Hence, precision measurements of particle physics experiments are expected to play important roles to search for the dark sector particles, as we will see in following sections.

Model parameters
There are eight parameters in the scalar potential: Since the two vacuum expectation values of the SM doublet H and the singlet Φ are fixed #5 The mediator φ seems to have unsuppressed interactions with h through the couplings A φH and λ φH even if the mixing angle is very close to zero. The coupling A φH is, however, suppressed when the mixing angle is small, as can be deduced from eq. (5). On the other hand, the coupling λ ΦH can be still large, though interactions between φ and SM fermions, which plays an important role in cosmology, are suppressed by small Yukawa couplings especially when the dark matter mass (hence, the freeze-out temperature) is small. and the Higgs mass is determined to be m h 125 GeV thanks to the LHC experiment, the number of free parameters in the potential becomes five. Adding the WIMP mass m χ as well as two couplings between the WIMP and the singlet Φ, namely c s and c p , the total number of free model parameters is then reduced to be eight in the end.
The condition of the two vacuum expectation values, v H 246 GeV and v Φ = 0, gives two relations among the model parameters, µ 2 H + λ H v 2 H /2 = 0 and µ 3 1 + A ΦH v 2 H /2 = 0, so that we can drop the two parameters µ 2 H and µ 3 1 from the set of independent parameters. Moreover, as seen in eq. (3), the Higgs mass is determined mainly by the parameter λ H when the mixing angle θ is suppressed, and hence it can also be dropped from the parameter set. On the other hand, when θ 1, the mediator mass m φ is determined mainly by the combination of the two parameters µ 2 Φ and λ ΦH , while the mixing angle is given mainly by the parameter A ΦH , as can be seen in eq. (3). Thus, we adopt m φ and sin θ as independent parameters instead of λ φH and A φH . As a result, we have the following eight parameters, m χ , c s , c p , m φ , sin θ , µ 2 φ , µ 3 and λ Φ , as free model input parameters.

Mediator decay
Because the decay of the mediator φ plays an important role in phenomenology of the model, we discuss some details of the decay in this subsection. As already mentioned in section 2.2, φ behaves like a light Higgs boson, and its decay into SM particles is from the mixing between H and Φ. Its partial decay width into a specific SM final state is We will discuss each decay mode below, and present how we estimate the decay width.
When the mediator is lighter than the electron-positron threshold, namely m φ ≤ 2m e , it decays mainly into two photons. We use the formula in Ref. [53] to calculate its partial decay width. Above the threshold, though this channel is not a dominant one anymore, the same formula is still used up to m φ = 0.6 GeV. If m φ ≥ 2 GeV, the decay width is computed by using the HDECAY code [54]. In the region of 0.6 GeV ≤ m φ ≤ 2 GeV, the formula in Ref. [55] is used to connect the regions m φ ≤ 0.6 GeV and m φ ≥ 2 GeV smoothly.
At the region of 2m e ≤ m φ ≤ 2m µ with m µ being the muon mass, it decays mainly into a electron and a positron. Its partial decay width is given by the following formula, When the mediator mass is above the muon threshold but below the pion threshold, namely 2m µ ≤ m φ ≤ 2m π , the mediator decays mainly into a muon pair. Its partial decay width is computed by the same formula as above with m e being replaced by m µ . The formulae for the two channels are used up to m φ = 2 GeV. In the mass region of m φ ≥ 2 GeV, the HDECAY code is used to compute the partial decay widths of the two channels.
When φ is heavier than the pion threshold but lighter than a few GeV, it decays mainly into a pair of pions (and a pair of K mesons if m φ ≥ 2m K ). Concerning the decay channel into ππ, we use the result in Ref. [56] to compute its partial decay width in the mass region of 2m π ≤ m φ ≤ 1.4 GeV, which is also consistent with the latest result in Ref. [57]. On the other hand, the width becomes negligibly small compared to other channels when m φ ≥ 2 GeV, so that we set Γ (φ → ππ) = 0 in this region. The width is evaluated by a linear interpolation in the mass range of 1.4 GeV ≤ m φ ≤ 2 GeV. On the other hand, the result in Ref. [56] is used again to compute the partial decay width of the φ → K K channel at 2m K ≤ m φ ≤ 1.4 GeV. Then, the width is evaluated using a linear interpolation at 1.4 GeV ≤ m φ ≤ 2 GeV by connecting the width to that of the φ → ss channel continuously, where s is the strange quark. The partial decay width of the φ → ss channel is computed by the HDECAY code at m φ ≥ 2 GeV. #6 Moreover, in order to take other hadronic decay channels into account, we also consider the φ → g g channel with g being the gluon in the mass range of m φ ≥ 1.4 GeV. When m φ ≥ 2 GeV, we use the result in the HDECAY code to compute its partial decay width, while the width is evaluated at 1.4 GeV ≤ m φ ≤ 2 GeV by a linear interpolation with the boundary condition of Γ (φ → g g) = 0 at m φ = 1.4 GeV.
Other decay channels open in the mass range of m φ ≥ 2 GeV (e.g. decays into a tau lepton pair, a charm quark pair, a bottom quark pair, etc.). We consider all possible decay channels in this region and compute their partial decay widths using the HDECAY code.
Here, we consider the φ decay into ππ (K K) in more details, because it is known to have a large theoretical uncertainty in the range of 2m π (2m K ) ≤ m φ ≤ 1.4 GeV [58] due to non-perturbative QCD effects. For instance, the total width of the SM Higgs boson in this mass range is predicted in other literature to be larger [59] or smaller [60] than what we have estimated based on Ref. [56]. In order to make our analysis conservative, we introduce a nuisance parameter σ to take this uncertainty into account. Then, the partial decay widths of φ → ππ and φ → K K channels are computed according to the equations where Γ ππ and Γ K K are the partial decay widths into ππ and K K computed based on Ref. [56].
In addition to the decay channels into various SM particles, the mediator particle φ can also decay into a pair of WIMPs when the mediator mass is larger than twice the WIMP mass, namely m φ ≥ 2m χ . Its partial decay width is described by All results obtained so far in this subsection are summarized in Fig. 1, where the total decay width of φ is shown in the left panel of the figure assuming that sin θ = 1 and φ does not decay into a pair of WIMPs. The gray band indicates the theoretical uncertainty due to non-perturbative QCD effects, which are taken into account by the nuisance parameter σ. In the right panel, the partial decay widths of several channels of φ are depicted. #6 In order to take the effect of the K meson threshold into account, we multiply the partial decay width of the φ → ss channel (computed in the HDECAY code) by the phase space factor of (m 2 A similar prescription is also applied for the other channel φ → cc with c being the charm quark. does not decay into a WIMP pair. The gray band indicates the theoretical uncertainty due to nonperturbative QCD effects. (Right panel) Partial decay widths contributing to the total width.

Higgs decay
In this model, the partial decay width of the Higgs boson into SM particles is slightly modified from the SM prediction because of the mixing between H and Φ. Its partial decay width into a SM final state (a fermion or a gauge boson pair) is given by the following formula: where the Higgs mass is fixed to be 125 GeV. In addition to these channels, the Higgs boson decays into a pair of WIMPs, and its partial decay width is given by the same formula as the one shown in eq. (12) with cos θ and m φ replaced by sin θ and m h , respectively: The Higgs boson can also decay into several mediators when the mediator φ is light enough. Considering the fact that the mixing angle θ is phenomenologically required to be much less than one and multi-φ channels are suppressed by their final state phase spaces, only the decay channel h → φφ can be potentially comparable to the other decay channels into SM particles (and a WIMP pair). The explicit form of its partial decay width is where the coefficient c φφh is given in eq. (7). This decay channel gives distinctive signals at the LHC experiment, depending on the mass and the decay length of the mediator. For instance, when the mediator decays mainly into leptons and its decay length is enough shorter than 1 mm, this Higgs decay channel gives a signal of h → 4 . On the other hand, if the decay length is much longer than the detector size, the channel gives a signal of the invisible h decay, regardless of the decay channel of the mediator. When the decay length is around the detector size, we can expect various signals at displaced vertex searches.  Table 1: Preselection criteria that we have imposed in our analysis. The second and third columns are for present and near future experiments/observations used to apply the criteria. The last column is for the section where each criterion is discussed in detail. The third column is now blanked, as the preselection criteria are not stronger than other constraints we discuss in following subsections.

Constraints
We introduce all constraints used in our comprehensive analysis. The likelihood of the constraints will be modeled in various functions which will be further discussed in section 4.1.
The usage and information of the constraints are discussed in the following in details.

Preselection criteria
Before performing the likelihood analysis, we apply preselection criteria on the model parameter space as one of our prior distribution. The following three criteria are imposed: CP conserving criterion, the criterion on the maximal value of the mixing angle, and the vacuum stability criterion. Those are implemented by a 1/0 logical cut in our analysis. All the preselection criteria involved in our analysis are summarized in Table 1.

CP conserving criterion
The pseudo-scalar interaction between the WIMP χ and the singlet Φ in eq. (1) breaks the CP symmetry under the general scalar potential of V (Φ, H). Then, the existence of the non-zero pseudo-scalar coupling c p induces the so-called s-wave WIMP annihilation into SM particles, namely the WIMP annihilates into SM particles without any velocity suppression at present universe. Moreover, when both the scalar and pseudo-scalar couplings, c s and c p , are nonzero, the WIMP annihilates into a pair of the mediators φ without the velocity suppression. In such cases, the WIMP annihilation cross section at present universe is almost the same as the one during the freeze-out process in the early universe, and this fact means that the cross section at present universe is expected to be about 1 pb.
WIMP annihilation cross section of (1) pb at present universe is, however, not favored by the CMB observation, when the WIMP mass is less than (1) GeV [61]. This is because the annihilation distorts the recombination history of the universe, while the observational result is consistent with the standard predication without the contribution from the WIMP. Thus, in order to avoid the CMB constraint, we set the pseudo-scalar coupling constant to be zero, namely c p = 0, in our analysis. In other words, we impose the CP symmetry on the interactions relevant to the WIMP. After switching off the coupling c p , the WIMP annihilates into SM particles as well as a pair of the mediators with a velocity suppression at present universe, which enables us to avoid the CMB constraint easily.

Criterion on |θ |
As already mentioned in previous section, the mixing angle θ is severely constrained by various experiments. In order to reduce a computational cost in our numerical analysis, we impose a condition on the angle as |θ | ≤ π/6, as a preselection criterion. This condition comes from the precision measurement of Higgs boson properties at the LHC experiment [62]. The measurement of the Higgs boson production process g g → h followed by its decay into a pair of gammas or W bosons gives a constraint on θ as cos θ ≥ 0.9, which is translated into the constraint on the angle θ as above. On the other hand, much severe constraints are eventually obtained from other particle physics experiments and cosmological observations, and those are taken into account in the subsequent likelihood analysis.

Vacuum stability criterion
The minimal WIMP model in eq. (1) has an extended scalar sector composed of Φ and H. In order to guarantee the stability of our electroweak vacuum, we impose the following condition on the scalar potential. We first define the fields ξ and η as Φ = η and H = (0, ξ/ 2) T , respectively, after taking the unitary gauge. Then, we impose a condition on the potential as V(η, ξ) ≥ V(v Φ , v H ) in the range of the fields |ξ| ≤ 1 TeV and |η| ≤ 1 TeV. Philosophy behind the condition is as follows: We consider the minimal model as an effective theory of the WIMP defined at the energy scale of 1 TeV, and require our electroweak vacuum to be absolutely stable within the range where the effective theory is applied. One might think that it is even possible to put a more conservative constraint on the potential by using the meta-stability condition, where our vacuum is required to be, at least, meta-stable with its lifetime much longer than the age of the universe. On the other hand, other constraints from particle physics experiments and cosmological observations give severer constraints on the potential, as we will see in the following subsections. #7 Hence, we adopt the absolute stable condition as a preselection criterion applied before the likelihood analysis.

Region of the parameters scanned
In addition to the the preselection criteria mentioned above, we scan the parameter space (m χ , c s , m φ , θ , µ 2 φ , µ 3 and λ Φ ) over the following ranges in our numerical analysis: The dimensionless coupling constants c s and λ Φ vary between −1 to +1, assuming that UV completion behind the minimal WIMP model is described by a weak interacting theory. On the other hand, dimensionful coupling constants basically vary within the energy scale of #7 Radiative corrections to the potential are also not taken into account because of the same reason. 1 TeV due to the validity of the minimal WIMP model. Upper limit on the WIMP mass m χ is fixed to be 30 GeV, simply because we are interested in the light WIMP region.
The surviving parameter space after applying the three preselection criteria are shown in appendix A. Though the seven parameters in eq. (16) are used to put the preselection criteria, #8 two of them, µ 2 Φ and λ Φ , are not very much relevant to the result of our likelihood analysis. This is because the coupling constant of the quadratic term, µ 2 Φ , appears only in the mass matrix of eq. (3), but the mass matrix is already parameterized by two parameters, m φ and θ . #9 The other parameter λ Φ appears in scalar quartic interactions, however no significant constraints on the quartic interactions are obtained so far and even in the near future. As a result, we will present our results of the likelihood analysis in terms of the following five parameters m χ , c s , m φ , θ and µ 3 in subsequent sections.

Conditions/Constraints from cosmology and astrophysics
In this subsection, we summarize cosmological and astrophysical conditions/constraints used to figure out the present status and future prospects of the minimal WIMP model. Those are taken into account through the likelihood analysis, unless otherwise stated. All the conditions/constraints involved in our analysis are summarized in Table 2.

Relic abundance condition
The WIMP should satisfy the so-called relic abundance condition. On observational side, the relic abundance of the WIMP, or in other words, the averaged mass density of the WIMP at present universe, is very precisely measured by the PLANCK collaboration [63]: with h being the normalized Hubble constant. The uncertainty of the observation is less than 2%, which is comparable to precise measurements at collider experiments.
On theoretical side, we calculate the relic abundance using the MicrOMEGAs code [77] based on the minimal WIMP model. The code first calculates the (thermal-averaged) annihilation cross section of the WIMP, and next compute the abundance by solving the Boltzmann equation numerically. Concerning the cross section, when m χ ≥ m φ , the WIMP annihilates #8 Five parameters (m φ , θ , µ 2 φ , µ 3 and λ φ ) among seven parameters are relevant to the preselection criteria. #9 Since µ 2 Φ determines the value of the parameter λ ΦH via eq. (3), it is scanned as a nuisance to be honest. mainly into a pair of the mediators, where the mediators eventually decay into SM particles. On the other hand, when m χ ≤ m φ , the WIMP annihilates into SM particles through the exchange of the mediator (or the Higgs boson) in the s-channel. Because the process is suppressed by small Yukawa couplings and the mixing angle, only the resonant region with m χ ∼ m φ /2 satisfies the relic abundance condition. Concerning the Boltzmann equation, it is known to have a theoretical uncertainty originating in the massless degrees of freedom of the universe during the freeze-out process. Since we are interested in the light WIMP scenario, the freeze-out temperature can be as low as the energy scale of the QCD phase transition, which makes the estimate of the massless degrees of freedom uncertain due to several non-perturbative QCD effects. Its uncertainty is reported to be at most 10% [78], and we adopt this maximum value to make our analysis conservative. Here, it is also worth mentioning that the mediator may contribute to the massless degrees of freedom and affect the Hubble expansion rate if it is light enough. However, we do not include this contribution in this analysis, because it is negligibly small compared to the above QCD uncertainty.
Note that it may have the Sommerfeld effect on the WIMP annihilation process [79][80][81], for the mediator is sometimes much lighter than the WIMP in our setup [82][83][84][85][86]. We found that the effect is not sizable for model parameter sets passing all conditions and constraints adopted in this paper. This is verified by computing the quantity, c 2 s /(4π) · (m χ /m φ ), and confirmed that it is always enough smaller than one for the sets. #10 We therefore do not take the effect into account to calculate the WIMP annihilation cross section.

Kinematical equilibrium condition
We also impose the kinematical equilibrium condition on the model, where the WIMP and SM particles are required to be in thermal equilibrium during the freeze-out process. Though the condition is automatically satisfied for a typical WIMP with its mass of (100) GeV, it should be imposed independently for the light WIMP, because both WIMP and mediator connect to SM with small couplings, and the condition is not automatically satisfied. It is worth emphasizing that we adopt this condition to figure out a very conventional WIMP parameter region in our setup. On the other hand, the condition can be relaxed by requiring that the WIMP is in the equilibrium at some temperature of the universe before the freezeout, because it still allows us to make a quantitative prediction on its abundance. We will discuss in appendix D how the result of our analysis alters by relaxing the condition.
The WIMP annihilates directly into SM particles when m χ ≤ m φ . Its reaction rate at the freeze-out temperature T f is estimated to be the product of the thermally averaged annihilation cross section (times the relative velocity) and the number density of the WIMP, namely Γ χχ ∼ 〈σ χχ v〉 T f n WIMP (T f ). The reaction rate Γ χχ becomes the same order as the expansion rate of the universe H(T f ) due to the relic abundance condition. The existence of the annihilation process guarantees that of the scattering process between the WIMP and the SM particles because of the crossing symmetry, and its reaction rate is Γ χSM ∼ 〈σ χSM v〉 T f n SM (T f ) with σ χSM and n SM (T f ) being the scattering cross section and the number density of the SM particles, respectively. Since the number density of the SM particles is much larger than that of the WIMP at the freeze-out temperature, this fact gives Γ χSM Γ χχ ∼ H(T f ) unless the annihilation cross section is significantly boosted compared to the scattering cross section. #10 This quantity is known to be the one for testing whether the Sommerfeld effect becomes sizable or not [80][81][82]. If it is smaller than one, the effect only gives a small correction to the thermally averaged cross section.
Hence, the equilibrium condition is usually automatically satisfied as far as m χ ≤ m φ . #11 On the other hand, the WIMP annihilates mainly into two mediators when m χ ≥ m φ , so that Γ χSM ≥ H(T f ) is not always guaranteed. However, even if Γ χSM ≤ H(T f ), the WIMP can be in the kinematical equilibrium with SM particles when the equilibrium is maintained between the WIMP and the mediator and between the mediator and the SM particles simultaneously. It means that, even if the reaction rate between the WIMP and SM particles is smaller than the expansion rate of the universe, the WIMP has a possibility to be in the kinematical equilibrium via the mediator in the universe. We thus take the following strategy to impose the kinematical equilibrium condition to the minimal WIMP model.
At each set of the input model parameters to define the model, we first calculate the freeze-out temperature T f by the MicrOMEGAs code. We next calculate the reaction rate Γ χSM and compare it with the expansion rate of the universe H(T f ). We accept the set if Γ χSM ≥ H(T f ). If it is not, we further calculate two reaction rates between the WIMP and the mediator and between the mediator and SM particles. The former reaction rate is estimated to be Γ χφ ∼ 〈σ χφ v〉 T f n φ (T f ) with σ χφ and n φ (T f ) being the scattering cross section and the number density of the mediator at the freeze-out temperature. The latter reaction rate has a more complicated form than the former one. In fact, three different processes contribute to the reaction rate; decay, scattering and absorption processes. The rate is estimated to σ φSM and σ φSM being the total decay width, the scattering and absorption cross sections of φ, respectively. If both Γ χφ and Γ φSM are larger than H(T f ), we accept the set of the model parameters. Note that, unlike the relic abundance condition, the kinematical equilibrium condition that we have discussed above is implemented by a 1/0 logical cut in our likelihood analysis.
In Fig. 2, the reaction rates Γ χχ , Γ χSM , Γ χφ and Γ φSM as well as the expansion rate of the universe H(T f ) are depicted as a function of the temperature of the universe in two cases; one is for the model parameter set that does not satisfy the kinematically equilibrium condition (left panel) and the other is for that satisfying the condition (right panel). The freeze-out temperature is shown as a vertical (orange) line in both panels. See also the figure caption for the model parameters used to calculate the reaction rates. As we can see from the comparison between the results in the two panels, the equilibrium condition gives the lower limit on the mixing angle |θ |. In order to manifest how the equilibrium condition works for the minimal WIMP model, the parameter region survived after imposing the equilibrium condition as well as the preselection criteria and the relic abundance condition is presented in appendix B. All explicit forms of the reaction rates are also given there.

Constraint from direct dark matter detections
Direct dark matter detection is known to be a stringent constraint for the WIMP based on the scattering between the WIMP and a nucleon. In the minimal WIMP model, the scattering occurs through the exchange of the mediator or the Higgs boson in the t-channel, and it contributes to the spin-independent scattering [85]. In our analysis, the scattering cross section is computed by using MicrOMEGAs code, where its explicit form is given by #11 When m χ is very close to m φ /2 and the decay width of the mediator is very suppressed, 〈σ χχ v〉 T f is indeed significantly boosted, and it requires a special treatment to calculate the correct relic abundance [87]. Here, m N is the mass of a nucleon and f N = f Tu .0191, f Ts 0.0447 and f T G 0.921, respectively. Since m φ m h is required to satisfy the relic abundance condition, the scattering process through the exchange of the mediator, namely the first term in the parenthesis of eq. (18), dominates the cross section.
On experimental side, the most stringent constraint on the spin-independent scattering cross section is from XENON1T [64], PANDAX [68,88], SuperCDMS [69], CRESST [65], Darkside-50 [72] and NEWS-G [71] experiments for the WIMP mass of our interest. On the other hand, in the near future, the constraint will be updated by LZ [73,74], Super-CDMS/SNOLAB [70] and NEWS-SNOLAB [66,67] experiments, if no signal is detected. In particular, the NEWS experiment utilizes several gas detectors (e.g. Helium to Xenon) and will play an important role for the search with the mass less than a few GeV. For the present constraints (XENON1T, CRESST, and Darkside-50), we use the Poisson distribution likelihoods given by the DDCalc code [43]. On the other hand, we involve the future constraints (SuperCDMS(SNOLAB), LZ, and NEWS-SNOLAB) assuming a half-Gaussian [89] form with the central value being set to be zero to figure out the future prospects of the light WIMP.

Constraint from N eff at T CMB
The mediator affects the expansion rate of the universe at the recombination era (T CMB 4 eV [75]) when it is lighter than the neutrino decoupling temperature, namely m φ T D (1) MeV [90]. This is because neutrinos are already decoupled from the thermal bath composed of photons and electrons when the temperature of the universe is below T D , while a part of the entropy of the universe is still being carried by the mediator at T ∼ T D and it is eventually injected into the two systems (neutrino and photon+electron systems) asymmetrically before the recombination according to the interaction of the mediator. #12 In the minimal WIMP model, the mediator interacts only with electrons and photons when its mass is small enough, so that the entropy carried by the mediator goes into those particles. As a result, this injection contributes to the photon temperature of the universe, which makes #12 If the mediator is lighter than T CMB , it contributes directly to the expansion rate of the universe at the recombination era and this possibility is already excluded by the constraint on N eff at the CMB observation. the difference between the photon and neutrino temperatures larger, namely the expansion rate of the universe at the recombination era becomes smaller than usual. #13 Such a contribution is severely constrained by the CMB observation, giving a lower limit on m φ in order not to alter the effective number of relativistic degrees of freedom N eff [75].
In the minimal WIMP model, the number N eff at the recombination era is predicted as [91] where [91]. Here, we have assumed that the mediator is never chemically and kinematically decoupled from the thermal bath, and it is indeed satisfied in most of the parameter region of the minimal WIMP model. This is because kinematical equilibrium is maintained among all species at the freeze-out temperature through the (inverse) decay process between the mediator and SM particles as well as the scattering process between the dark matter and the mediator, unless the mixing angle is vastly suppressed. #14 As a result, the mediator is in the equilibrium with SM particles below the freeze-out temperature.
On the other hand, in other parameter regions, the mediator is chemically decoupled from (but still maintains the kinematical equilibrium with) the thermal bath before the freeze-out of the WIMP and/or kinematically decoupled from (and eventually recoupled with) the thermal bath after the freeze-out. In these cases, the mediator can contribute to N eff more than that discussed in eq. (19) and receives a severer constraint from the observation. We, however, keep using the formula (19) as a conservative constraint obtained from the CMB observation, and leave those precise computations for a future study.
The effective number of relativistic degrees of freedom is observed to be N eff = 2.99 ± 0.17 by the Planck collaboration [92]. The data could be improved by future observations, which resolves small angular scales (high multipole number ) of CMB and allows a more precise measurement of N eff by observing its damping tail at high-. According to their forecast, the measurement of N eff can be improved by one order of magnitude at, for example, the CMB-S4 experiment [76], which leads to the expected limit of ∆N eff < 0.017, if we do not see any deviation from the standard prediction. We involve these (expected) limits in our analysis for the present status and future prospects of the minimal WIMP model.

Constraint from BBN
If the mediator decays during or after the era of the Big Band Nuclear synthesis (BBN), it may spoil this successful scenario of generating light elements in the early universe. Hence, the mediator should be constrained so that it does not spoil the BBN scenario, which depends on its mass, lifetime, decay channels and abundance at the time that it is decaying.
For instance, if the mediator is heavier than 2m π and mainly decays hadronically, it modifies the neutron-to-proton ratio and increases the primordial helium mass fraction, which gives a constraint on its lifetime as τ φ (1) second [93]. On the other hand, the mediator also affects the BBN scenario even when its mass is in the range of 4 MeV ≤ m φ ≤ 2m π #13 One may think that the WIMP also affects the expansion rate if it is light. The WIMP mass is, however, constrained to be more than 10 MeV due to the kinematical equilibrium condition and collider constraints as we will see later, so that the entropy carried by the WIMP is injected into the thermal bath at much above T D . #14 For instance, the mixing angle of (10 −5 ) is enough to maintain the kinematical equilibrium even if the mediator is lighter than (100) MeV, and it is not excluded by collider constraints as seen in next section. and thus decays leptonically. This is because a few MeV-electrons pass their energies to CMB photons through the inverse Compton scattering and generate MeV gammas, which disintegrate 2 H and 7 Be (which eventually forms 7 Li). As a result, the leptonic decay produces less 7 Li and more 3 He/ 2 H than the usual case, leading to a constraint as τ φ (10 5 ) seconds [93]. #15 The mediator may affect the BBN scenario by other mechanisms if it is very light. For instance, when m φ T BBN ∼ 0.1 MeV, the mediator is relativistic during the BBN era and contributes to the expansion rate of the universe as an additional light degree of freedom. On the other hand, even when the mediator is in the mass range of T BBN m φ T D , the mediator alters the expansion rate of the universe, because it changes the ratio between the photon and neutrino temperatures as discussed in section 3.2.4.
In order to impose the BBN constraint on the minimal WIMP model in a rigorous way, we have to consider the cosmology of the mediator at the late universe [94,95] #16 and implement all the decay channels of the mediator in an appropriate BBN code, which is beyond the scope of this paper. Instead, we impose the following two constraints, τ φ ≤ 1 second when m φ ≥ 2m π and τ φ ≤ 10 5 seconds when 4 MeV ≤ m φ ≤ 2m π , by a simple 1/0 logical cut in the likelihood analysis, which is indeed a reference constraint adopted in many literature (see, e.g. Ref. [52]). #17 We leave precise computations of the BBN constraints for a future study. On the other hand, we do not involve the BBN constraints concerning the expansion rate of the universe, for those in section 3.2.4, ∆N eff had given severer constraint.
Let us briefly discuss how the BBN constraint works for the minimal WIMP model. The lifetime of the mediator is proportional to the mixing angle squared, and thus it gives the lower limit on the angle sin θ . As can be seen in Fig. 8 of appendix B.2, the result on the (m φ , sin θ )-plane, no lower limit on the mixing angle is obtained by imposing preselection criteria as well as relic abundance and kinematical equilibrium conditions, so that the BBN constraint plays an important role to make the allowed model parameter space finite.

Other constraints
It is also possible to impose further cosmological and astrophysical constraints on the minimal WIMP model in addition to those discussed above. One of such constraints is from the distortion of the CMB spectrum due to the late time decay of the mediator. When the lifetime of the mediator is longer than 10 6 seconds, where the double Compton scattering process (γ e → γ γ e) is not active and the CMB cannot maintain the Planck distribution against an additional photon injection, the constraint on the so-called µ-distortion parameter put an upper limit on the lifetime [96]. When the lifetime is longer than 10 9 seconds, where the Compton process (γ e → γ e) is also inactive and even the Bose-Einstein distribution cannot be maintained against the injection, the constraint on the y-distortion parameter put the upper limit on the lifetime [97]. These limits are, however, weaker than CMB and BBN constraints in Sec. 3.2.4 and 3.2.5. Hence, we do not include the constraints from the distortion of the CMB spectrum due to the late time decay of the mediator in our analysis.
The other possible constraint can be obtained by the observation of neutrinos from the supernova (SN) 1987A [98]. Physics behind the constraint is as follows. The collapse at #15 When the mediator is lighter than 4 MeV, the energy of the electron produced by the mediator decay is not sufficient to photo-disintegrate nuclei, so that the constraint discussed here cannot be applied [94]. #16 In the references, sin θ is assumed to be so small that the absorption and decay processes do not work to maintain the chemical equilibrium. Hence, constraints derived there cannot be directly applied to our case. #17 We would also like to thank Dr. Sebastian Wild for fruitful discussion for this BBN constraint [94].
the SN heats up the core and its energy is estimated to be (10 59 ) MeV, so that the cooling rate must be, at least, smaller than this value. The neutrino emission during a supernova explosion is known to be the only mechanism to cool down the core within the SM. On the other hand, the mediator can also contribute to the cooling in the minimal WIMP model. The dominant process is the nucleon scattering, N N → N N φ, followed by the φ decay into SM particles. It is then possible to put an upper limit on the mixing angle sin θ , such that the instantaneous luminosity in novel particles does not exceed the value (10 52 ) erg/s. When sin θ is too large, the parameter space is allowed because φ decays or is trapped inside the core and not contribute to the cooling. As a result, the SN physics has a potential to put a constraint on a certain region of m φ and sin θ .
To include SN1987a cooling constraints into our analysis, it can be a highly non-trivial task because SN is a complicated physics system and the constraint from it suffers from several uncertainties [99][100][101][102][103][104]. For example, the nature of the core of protoneutron star to the primary driver of the shock revival, the temperature, the density profiles, and equation of state of the progenitor star and the cross section of various QCD processes with soft radiations and environment effects to the process are known to give such uncertainties. Unfortunately, none of studies concerning the constraint is found to include all the uncertainties, and such a comprehensive study of the constraint is beyond the scope of our study. Hence, we did not include the constraint in our likelihood analysis, while we present the SN1987A exclusion contour on top of our result in the (m φ , sin θ )-plane in section 4.2.
Finally, we would like to comment on the cross section of the dark matter self-scattering process (χχ → χχ). Since the mediator φ in this model can be a thousand lighter than the dark matter and leads to an interesting enhancement of the dark matter self-scattering cross section. This enhancement becomes significant and velocity-dependent when the incident dark matter velocity is as small as 10 −3 because of the diagram exchanging a light mediator in the t-channel. On the other hand, the velocity-dependence disappears when the velocity is larger than (10 −2 ) and the self-scattering cross section is approximately given as [83] Phenomenologically, the velocity-dependent self-scattering at v ∼ (10 −3 ) could provide a solution to the "small scale crisis" [105], namely the so-called core-cusp problem of dwarf galaxies can be solved by the self-scattering. On the other hand, a robust constraint on the self-scattering with v ∼ 10 −2 is obtained from various clusters of galaxies (see, eg., Ref. [106]). We have checked that all parameter points satisfying all conditions and constraints imposed in this paper are also consistent with the constraint on the self-scattering, and thus we do not include this self-scattering constraint in our likelihood analysis.

Constraints from collider experiments
In this subsection, we summarize collider constraints that we have considered for our likelihood analysis. Since the WIMP couples to SM particles mainly through the mediator while the mediator couples to several SM particles, the constraints on the minimal WIMP model are mostly from mediator productions at various colliders. #18 Collider signals then depend #18 Only exception comes from the interaction between the WIMP and the Higgs boson. This interaction is constrained by observing the invisible decay width of the Higgs boson at colliders, as seen in section 3.3.4.

Present
Future Section Υ decay CLEO [107], BABAR [108,109] Belle II [110] 3.3.1 B decay Belle [111,112], LHCb [113][114][115]  on how the mediator is produced and how it decays, depending on its mass m φ and the mixing angle sin θ . We will discuss below the constraints based on the production processes.
All the collider constraints that we have involved in our analysis are summarized in Table 3.

Upsilon decay
Upsilons can produce the mediator through their decays. The most important decay channels are Υ (1S, 2S, 3S) → γ φ followed by the φ decay into a lepton pair such as µ − µ + and τ − τ + . These upsilon decays allow us to detect the mediator through a narrow peak search in the di-lepton invariant mass spectrum at the mass region of m φ < m Υ 10 GeV. These channels have been searched for by CLEO [107] and BaBar [108,109] collaborations. Since no signal has been detected so far, they put an upper limit on the following quantity: G F , m b and α are the Fermi constant, the bottom quark mass and the fine structure constant, respectively. The parameter QCD takes a value between 0.5 to 1.5, originated from QCD bound state and relativistic corrections. The branching fraction of the decay channel Υ → µ − µ + have been measured precisely [140]. Then, an upper limit on Br(Υ → γφ) × Br(φ → ) is obtained at 90% confidence level, depending on the mediator mass m φ [107][108][109]. For instance, an upper limit on Br(Υ → γφ) × Br(φ → µ − µ + ) 3 × 10 −6 is obtained in the mass region of m φ 3.5 GeV, while Br(Υ → γφ) × Br(φ → τ − τ + ) 3 × 10 −5 is obtained for 4 GeV m φ 8.5 GeV. The limit becomes weaker when m φ approaches to m Υ .
The mediator is assumed to promptly decay into a lepton pair in the above analysis. One might suspect that this method does not work because the mixing angle may become so small that the mediator is too long-lived. However, the method indeed works well but its reason is a bit complicated. We take the mediator mass to be m φ = 220 MeV as an example, just above the µ − µ + threshold giving the longest lifetime to make φ decay into µ − µ + at each mixing angle. #19 The decay width of the mediator is then estimated to be Γ φ 10 −9 × sin 2 θ GeV. The constraint Br(Υ → γφ) × Br(φ → µ − µ + ) 3 × 10 −6 is then translated to sin θ 0.14, so that the decay length of the mediator produced by the Υ decay is estimated to be γ φ τ φ c = (0.1) mm with the speed of light c. Here, the so-called gamma factor is estimated to be γ φ m Υ /(2m φ ) 25. This decay length is much shorter than the present sensitivity of displaced vertex searches, which requires (1) cm [117]. #19 Needless to say , the lifetime (thus, the decay length) becomes shorter if the mediator mass is heavier. The Belle II experiment will update the constraints in the near future with the integrated luminosity about fifty times more than that of the BaBar experiment, if no signal is detected there [110]. The constraints on the branching fractions will be about six times more severer than the present ones. As a result, for instance, the constraints will be Br(Υ → γφ)×Br(φ → µ − µ + ) 5 × 10 −7 and Br(Υ → γφ) × Br(φ → τ − τ + ) 5 × 10 −6 . Note that the decay length of the mediator is, at most, (1) mm even in the near future, so that its decay can be treated as a prompt one. We include all the constraints discussed here in our likelihood. It, however, turns out that the present constraints are not stronger than others. Those from B(K) decays are more sensitive in the mass region of m φ < m B (m K ) with m B (m K ) being the B(K) meson mass, while the LEP experiment gives a stronger constraint in m B < m φ < 9.2 GeV, as we will see in following discussions. On the other hand, the future constraints could play an important role to search for the mediator in a certain mass region of φ.

B meson decay
When the mediator φ is lighter than B mesons (m B 5.3 GeV), it can be produced by their decays through the sub-process b → sφ, where it is induced mainly from the one-loop diagram composed of weak bosons and various quarks in the loop [52]. Among various B meson decays, one of the efficient channels to search for the mediator at collider experiments is the charged B meson decay, B ± → K ± φ. Its decay width is estimated to be [141]: where the scalar form factor of K ± is estimated to be F K (q) 0.33 (1 − q 2 /38 GeV 2 ) −1 , while the coefficient of the FCNC effective interaction, C sbsL b R φ + h.c., is given by |C sb | |2g 2 m b m 2 t V * ts V t b sin θ |/(64π 2 m 2 W v H ) 6.4×10 −6 sin θ with V t b and V ts being the CKM matrix elements. Combined with the total decay width of the B meson, Γ B + 4.1 × 10 −13 GeV, the above formula allows us to compute the branching fraction of the decay channel B ± → K ± φ. It is also possible to consider other decay channels for the φ search such as B 0 → K * 0 φ. Collider signals then depend on how the mediator decays and how long distance it travels before the decay. In this subsection, all experimental results of B meson decay at present and expected null signal results in the near future are discussed.

B meson decay with prompt φ decay
When the decay length of the mediator is enough smaller than the detector size, its decay can be treated as a prompt one. In such a case, the decay channel B ± → K ± φ → K ± µ − µ + is frequently used to detect the mediator, as it gives a clean signal against backgrounds. At present, BaBar [118], Belle [111] and LHCb [113] collaborations give upper limits on its branching fraction as P p Br(B ± → K ± φ) Br(φ → µ − µ + ) 3 × 10 −7 , accepting that the theoretical prediction of the SM contribution is Br(B ± → K ± µ − µ + ) SM = (3.5 ± 1.2) × 10 −7 [142].
Here, the prefactor P p is the probability that the mediator decays inside the detector, where τ φ is the lifetime of the mediator, while θ φ is the angle between the direction of the mediator produced from the B decay and that of the beampipe. The boost factor is estimated to be γβ m B /(2m φ ), as the mediator is produced from the B decay. The size of the detector is set to be l x y 25 cm referencing to those of Belle and BaBar detectors [58].
The theoretical error on the SM contribution overwhelms the experimental errors due to the uncertainty of the form factor in eq. (22). Therefore, a significant improvement of the sensitivity is not expected even if more data is accumulated at e.g. the upcoming SuperKEKB experiment [116]. On the other hand, if the theoretical error is much reduced, thanks to the development of lattice QCD simulations and/or the use of the characteristic feature coming from the narrow dilepton invariant mass peak from the mediator decay, the analysis can have a further improvement. Since no convincing study on the issues has been performed yet, we do not involve the constraint from the prompt B ± → K ± φ → K ± µ − µ + decay to investigate the future prospects of the minimal WIMP model and involve it only for its present status.

B meson decay with displaced φ decay
When the decay length of the mediator φ becomes comparable to the detector size, displaced vertex is the most powerful channel to search for φ. At present, BaBar [117] and LHCb [114,115] collaborations give various stringent constraints. The BaBar collaboration has looked for the displaced vertex caused by decay channels φ → e − e + , µ − µ + , π − π + and K − K + by observing the invariant mass spectrum of the decay products, where the mediator is produced mainly from B decays, B → X s φ with X s being a hadronic system with a strangeness. Since no excess over the SM prediction has been observed, they put a constraint on Br(B → X s φ) Br(φ → e − e + , µ − µ + , π − π + , K − K + ), which is translated to that on the mixing angle, sin 2 θ 2 × 10 −8 at 90% C.L., almost regardless of m φ and τ φ in the region of 0.5 GeV ≤ m φ ≤ 1.5 GeV and 1 cm ≤ cτ φ ≤ 20 cm. We adopt it in our likelihood analysis. More stringent limits on the mixing angle for regions m φ ≥ 1.5 GeV and cτ φ ≥ 20 cm is obtained at LHCb (see below) and beam dump experiments, respectively.
On the other hand, the LHCb collaboration is recently searching for the mediator produced by both charged and neutral B-meson decays: B ± → K ± φ → K ± µ − µ + and B 0 → K * 0 φ → K * 0 µ − µ + . Their latest report based on the negative result of the search using 3 fb −1 data at 7 TeV and 8 TeV running is found in Ref. [115]. We also involve this result in our analysis, where it gives a stringent constraint in the region of m φ ≥ 1.5 GeV.
In the near future, Belle II collaboration will update the constraint (given by the BaBar collaboration) with 100 times more data (50 ab −1 ) [116], if no signal is detected. On the other hand, the LHCb collaboration will accumulate data during the 13 TeV running, where 300 times more B mesons will be in data compared to the present one [121]. For the near future prospects in our analysis, we therefore adopt the BaBar and LHCb constraints mentioned above with their sensitivities on the mediator searches increased 10 and 17 times.

B meson decay with very long-lived φ decay
When the decay length of φ is much longer than the detector size, it is searched by utilizing the channel, B ± → K ± + missing, where the SM process, B ± → K ± νν, is a background against the signal. This channel is also used to search for the mediator decaying invisibly, B ± → K ± φ → K ± χχ, even if it is short-lived. At present, Belle [112] and BaBar [119,120] collaborations put a constraint as P l Br(B ± → K ± φ) + P p Br(B ± → K ± φ) Br(φ → χχ) ≤ 1.6 × 10 −5 . The factor P l is the probability that φ decays outside the detector, where the values of l x y and γβ are the same as those in eq. (23). In future, the Belle II collaboration will update the sensitivity as P l Br( . We include this constraint in our analysis for the near future prospects.

Kaon decay
The mediator is also produced from Kaon decays if it is lighter than Kaons (m K 0.5 GeV), which is induced from the sub-process s → d φ. The most efficient channels to search for φ are K ± → π ± + φ and K L → π 0 + φ, where their decay widths are given as follows [52]: Since scalar form factors for pions are close to unity [143], we neglect those in the above equations. On the other hand, the coefficient C ds comes from the FCNC effective interaction, 2+0.5i)×10 −9 sin θ with the CKM matrix elements V ts and V t d . By given the total decay widths of the Kaons, Γ K ± = 5.3×10 −17 GeV and Γ K L = 1.286×10 −17 GeV, we can compute their branching fractions. The Kaon decays are then followed by the φ decay, and the decay channels φ → µ − µ − , e − e + are used, as it gives a clean signal and has less theoretical uncertainties than others.
In the near future, an improvement of the sensitivity on the search with the K ± decay is not expected, for the systematic error at the scalar form factor already dominates. Moreover, there is no successor of the KTeV experiment, so that the constraints from the neutral Kaon decays are not expected to have any improvement in the near future. We hence do not consider the constraints relevant to the above prompt Kaon decays in our likelihood analysis for the near future prospects of the minimal WIMP model.

Kaon decay with displaced φ decay
When φ decays with the decay length of (100) m, it is efficiently searched at proton beam dump experiments. At present, the CHARM collaboration gives the stringent constraint on this search [128], where many Kaons and B mesons are produced from the 400 GeV proton beam which is dumped into the copper target. Then, the mediator is expected to be promptly produced from the channels, K ± → π ± φ, K L → π 0 φ and B → X s φ. At the experiment, the detector is located at 480 m away from the target with its size of 35 m, and φ penetrating the wall is searched for with leptonic decay channels φ → e − e + and µ − µ + .
Since zero signal event was observed as expected by background, the collaboration put a constraint on the number of signal events as N dec ≤ 2.3 at 90% C.L. [129]. The total number of signal events at the experiment is estimated by [57] where N p.o.t is the number of protons on the target, and n K(B) is the number of K(B) meson created per each incoming proton. Since Kaons are long-lived and absorbed in the target, the factor H /(cγ K τ K ) is multiplied to its contribution with γ K , τ K and H being the Lorentz factor, lifetime of the K meson and the hadronic absorption length, respectively. The factors P K dec and P B dec are probabilities that the mediator decays inside the detection region, with η K(B) geom and η K(B) rec are the geometric and reconstruction efficiencies, respectively, assuming those are independent along the beam line. In the CHARM experiment, the values of the above parameters are In the near future, the above constraint will be improved by the SHiP experiment at CERN SPS [131] if no signal is detected, where the 400 GeV proton beam is dumped to the fixed target and the expected number of protons on the target is N p.o.t = 2 × 10 20 . The SHiP detector is located at 69 m away from the target with its size of 51 m, so that the probability that the mediator decays inside the detection region is given by P dec in eq. (28) with η K  7). Then, the constraint will be updated as N dec ≤ 3 at 95% C.L., if no signal is detected there. Since the SHiP experiment can detect the decay channels of the mediator into ππ and K K, it will have more sensitivity for a heavy mediator ( 4 GeV) than that of the CHARM experiment. We include this future-expected limit to study the future prospects of the minimal WIMP model.

Kaon decay with very-long lived φ decay
When the lifetime of the mediator is very long, it is searched through the process, K ± → π ± +missing, where the SM process, K ± → π ± νν, becomes a background against the signal. In addition to the signal channel, K ± → π ± (long − lived φ), another channel, K ± → π ± φ → π ± χχ, also contributes to the signal when the mediator is twice heavier than the WIMP, as in the case of the B meson decay described at the last paragraph of section 3.3.2. At present, the E949 collaboration put a stringent constraint on the branching fraction P l Br(K + → π + φ) + (1 − P l ) Br(K + → π + φ) Br(φ → χχ), where P l is defined in eq. (24) with the size of the detector and the boost factor being replaced by l x y 145 cm and γβ 1, respectively. The constraint on the branching fraction at 90% C.L. is found in Fig. 18 of Ref. [125].
On the other hand, such a very long-lived mediator is also searched by using the neutral Kaon decay, K L → π 0 + missing. Then, the SM process, K L → π 0 νν, becomes a background against the signal. At present, the KEK E391a experiment put a constraint on the branching fraction as P l Br(K L → π 0 φ) + (1 − P l ) Br(K L → π 0 φ) Br(φ → χχ) ≤ 2.6 × 10 −8 with l x y and γβ m φ used in P l being replaced by l x y 1 m and γβ m φ 1 GeV, respectively [130]. We involve the above two constraints from the charged and neutral Kaon decay experiments in our likelihood analysis to investigate the present status of the minimal WIMP model.
In the near future, the constraint from the charged Kaon decay will be improved by the NA62 experiment as P l Br(K + → π + φ) + (1 − P l ) Br(K + → π + φ) Br(φ → χχ) 10 −11 if no signal is detected, where P l = exp[−l z /(γβ cτ φ )] with l z = 65 cm and γβ m φ = 37.5 GeV [126,127]. On the other hand, the constraint from the natural Kaon decay will be improved as P l Br(K L → π 0 φ) + (1 − P l ) Br(K L → π 0 φ) Br(φ → χχ) 1.46 × 10 −9 at the KOTO experiment with P l being the same as that for the KEK E391a experiment [132], which is close to the so-called Grossman-Nir bound. Since the branching fractions of the decay channels, Br(K + → π + φ) and Br(K L → π 0 φ), are at the same order, the constraint from the charged Kaon decay is stronger than the neutral one in the most of the parameter region. Despite this fact, we involve both the constraints in our likelihood analysis for the sake of comprehensiveness to investigate the future prospects of the minimal WIMP model.

Higgs decay
The mediator can also be produced from the Higgs decay, h → φφ, and this process is indeed being investigated by measuring the Higgs boson property carefully at the LHC experiment. The latest result of the measurement is consistent with the SM prediction, so that we have various constraints on the decay process depending on the decay length of the mediator.

Higgs decay with prompt φ decay
When the mediator promptly decays into two leptons, it is searched for through the Higgs decay channel, h → 4 . At present, the ATLAS collaboration put a constraint on the branching fraction of the four muon process, (P p ) 2 Br(h → φφ) × Br(φ → µµ) 2 , where P p is the probability that the decay is considered to be a prompt one at the experiment and defined as P p = 1 − exp[−σ/(γβ cτ φ )] with σ = 1 mm and γβ = m h /(2m φ ), respectively. The limit on the branching fraction at 90% C.L. is found in Ref. [134]. On the other hand, the CMS collaboration put constraints on similar branching fractions with the same σ and γβ by considering various leptonic channels, h → 4µ, 2µ2τ and 4τ. Constraint on the branching fraction of each channel is obtained from Ref. [133] using the same prefactor P p . We involve these constraints in our analysis for the present status of the minimal WIMP model. #20 In the near future, both the collaborations will update their constraints at the highluminosity upgrade of the LHC experiment (HL-LHC) [137], if no signal is detected. Since we can expect about a hundred times more data at the HL-LHC than the current one, the sensitivity of the search will be improved by one order of magnitude if the statistical error dominates. For the future prospects of the model, we hence involve the same constraints on the above branching fractions in our analysis with making the limits ten times severer.

Higgs decay with displaced φ decay
When the mediator decays with the decay length of 0.1-10 m, it is searched for by the displaced vertex analysis at the LHC experiment. The mediator is produced mainly from the Higgs decay, pp → h + X → φφ + X , when the mixing angle, sin θ , is small. Then, two displaced vertices are formed by the mediator decay if it is long-lived, and those are detected via the decays φ → µµ, ee and ππ. Here, the decay products are collimated as the mediator is boosted. At present, the ATLAS collaboration put a constraint on the branching fraction of the Higgs decay with displaced vertices using various decay channels of the mediator [144]. These results are summarized in Ref. [135], where an upper limit on the branching fraction, Br(h → φφ), is given at 90% C.L. as a function of the mediator mass and the mixing angle (lifetime). Roughly speaking, the region of Br(h → φφ) 30% is excluded when 0.3 GeV m φ 60 GeV and 0.1 m cτ φ 10 m. We involve the constraint in Ref. [135] in our analysis for the present status of the minimal WIMP model. Note that this constraint is not applied in the m φ ≥ 2m χ region, as the φ → χχ decay was not considered in Ref. [135].
Since the systematic error already dominates the statistical one to put the constraint on the branching fraction Br(h → φφ) (with displaced vertices) at the present LHC experiment, it is difficult to expect a significant improvement of the sensitivity on this search at the HL-LHC experiment. We therefore do not consider the constraint from the search for the Higgs decay with displaced vertices to study the future prospects of the minimal WIMP model.

Higgs decay with very long-lived φ decay
When the mediator is very long-lived, the decay process, h → φφ, contributes to the invisible decay width of the Higgs boson. In addition, the other decay channels exist, which are h → χχ and h → φφ → 4χ, and always contribute to the invisible decay width of the Higgs boson without respect to the decay length of the mediator. As a result, the new physics contribution to the branching fraction of the invisible Higgs decay is given as follows: 30 are the probabilities that the mediator decays outside and inside the detector, respectively, with the size of detector and the boost factor being 30 = 30 m and γβ m h /(2m φ ). At present, the constraint on the branching fraction is given as Br(h → inv.) BSM ≤ 0.19 at 90% C.L. [136] from Higgs precision measurements, and we include this in our analysis for the present status of the minimal WIMP model. #20 Both ATLAS and CMS collaborations assumed in their analyses that the Higgs boson decays into two pseudo-scalars, which means that different branching fractions of Higgs and mediator decays are used. We have corrected those to apply their constraints to our case where the Higgs boson decays into two scalars.
The constraint on this branching fraction will be improved at the HL-LHC experiment with 3 ab −1 data. If no signal is detected there [138], we can set Br(h → inv.) BSM ≤ 0.05 at 90% C.L. We include it in our analysis for the future prospects of the minimal WIMP model.

Direct production
When the mediator is heavier than the B meson and the mixing angle is not very suppressed, the mediator is also efficiently searched for by its direct production at high-energy colliders. The mediator is produced by the Bremsstrahlung process, ff → V * → V φ, with f and V being a SM fermion and a weak gauge boson, respectively, and the mediator always decays promptly in the parameter region of our interest. At present, the LEP experiment put the most severe limit on the mixing angle (sin θ ) as a function of the mediator mass (m φ ). To be more precise, the L3 collaboration have used the hadronic decay channel of the mediator associated with Z boson decays, Z → νν, e − e + and µ − µ + , where its result is found in Ref. [139]. Similarly, the ALEPH collaboration have also used the hadronic decay channel of the mediator, but only associated with the decay Z → νν. Its result is found in Ref. [145], which is slightly weaker than that of the L3 collaboration. On the other hand, the OPAL collaboration have utilized the inclusive Z production with its leptonic decays, Z → e − e + and µ − µ + , where its result is found in Ref. [146]. It is, however, applicable for a very light mediator, say as light as 10 −6 GeV, and it dose not play an important role in the analysis.
The LHC experiment is also searching for the mediator produced by the above direct production process [147], though it is less sensitive than the LEP experiment. Moreover, the obtained limit on the mixing angle is already dominated by systematic errors, it seems difficult to expect a significant improvement even at the HL-LHC experiment. We therefore involve the limit obtained by the L3 collaboration in our likelihood analysis to investigate the present status of the minimal WIMP model, and do not consider the expected constraint from the direct mediator production for the near future prospects of the WIMP model.

Results
We are now in a position to present the result of our likelihood analysis. We first explain the framework of our analysis in some details, as it plays a crucial role in our study. Then, we will present various results obtained in our analysis and discuss their implications to investigate the present status and near-future prospects of the minimal WIMP model.

Simulation framework
There are many released limits and data used in this work. Because of lacking discovery of new physics, the majority of them are based on the 90 % upper or lower confidence limit on a one-dimensional physical observable. #21 Nevertheless, for some other constraints such as relic density, because of its clear discovery, it is nature to describe such a likelihood probability by a two-tail Gaussian with a narrow peak. Except above two types of limits, theoretical constraint is used to veto the parameter space whose answer is often physical or not physical. As summarized in Table 4, we have used four types of likelihood function in our #21 The terminology 90 % C.L. for δχ 2 = 2.71 in one-tail Gaussian likelihood distribution is sometimes used.  Step function likelihood. In the following, we will present the usage of these four likelihood functions.
The most powerful advantage of the global analysis is to increase the statistics by combining different data sets which allows us to remove more parameter spaces from different corners. However, in order to conservatively exclude parameter space by adding the statistics, a proper likelihood function is needed and its tail at the exclusion region is particularly important. Once the central value and error bar are given, the Gaussian likelihood is On the other hand, for counting experiments, if expected event number and observed event number are provided, Poisson distribution is usually adopted for its likelihood function, where s and b are expected signal and background events, while o is observed events.
It is usually having null signal detection in new physics search. In such a search, a lower limit or upper limit of an observable is reported. Such a limit with some information given by experimental collaboration, one can of course reconstruct and verify some likelihood used in the experimental collaboration. Because of expansive background simulation and likelihood computation time, it will not be realistic to repeat whole procedure to reconstruct such a precise likelihood as used in the experimental collaboration. Therefore, we model the likelihood based on two reasonable assumptions: null signal detection and Gaussian distribution with the alignment of 90 % C.L. The likelihood of such a distribution called half Gaussian is almost identical as eq. (30) but with two differences; The center value is set to be zero because of null signal detection and the error-bar-squared can be obtained by the 90 % C.L. limit divided by 2.71 to align 90 % Gaussian one tail limit. Note that such a half Gaussian distribution shall be more conservative than precise one because new physics signal is expected to have an excess than background in the future.
Some experimental upper/lower limits are built based on multi-dimensional variables which are hard to reconstruct their likelihoods. For example, the limits of the LHC Higgs displaced vertex search in Sec. 3.3.4 are based on three dimensional variables, m φ , sin θ and the branching fraction of h → φφ. For a sake of simplicity, one may just use a hard cut for such an experimental constraint to answer whether the parameter space is allowed or excluded. Sometimes, one could perform a scan with the constraints described in case (ii) by using a step function, but it is particularly avoided in our analysis because it will lose the advantage of combined analysis as discussed in the previous paragraph.
Finally, when performing two-dimensional contour figures in this paper, we use the method "Profiled Likelihood", namely the minimum Chi-squared method. Like frequentist, such a method allows us to get rid of unwanted/not interesting parameters by taking the maximum likelihood along the direction of other unwanted/not interesting dimensions.
Regarding to the sampling, as seen in Appendices A and B, the likelihood function is flat within vast 1sigma allowed region but dramatically changing at the 2sigma tails. Therefore, in our analysis, we adopt an unusual strategy to scan our parameter space. We first scan the parameter space with the range in eq. (16) and the conditions in Table 1. Note that two dark matter parameters m χ and c s are not scanned in this step. After collecting (10 6 ) allowed points, we use it to present the parameter space constrained by the preselection criteria in Appendix A. In the second step, based on the former collected data set, we varied m χ and c s for each points to fit the relic abundance and kinematical equilibrium conditions. However, the parameter λ Φ does not affect the dark matter phenomenology at all. In this step, we evaluated (10 8 ) points but only ∼ 5×10 6 points fulfill the relic abundance and kinematical equilibrium conditions. In the third step, based on the (10 8 ) points collected from the former scan, we evaluate all the likelihoods with the constraints given in Tables 2 and 3. Finally, we use the knowledge we learned from previous scans as the new prior distribution to do several more sophisticated scans. We employed two scan tools (emcee [148] and MultiNest [149]) to undertake the sampling. With total (10 8 ) likelihood evaluation done by the tools emcee and MultiNest, we found the coverage is good enough. We combine all the scans from the former steps to draw our figures.

Present status
Results of our analysis for the present status of the minimal WIMP model projected on the (m χ , m φ )and (m φ , | sin θ |)-planes are shown in the left and right panels of Fig. 3, respectively, while those projected on all planes of input parameters are found in appendix C.1.

(a)
In the left panel, the upper limit on the mediator mass m φ comes from the relic abundance condition. As mentioned above, most of parameter regions satisfy the relation m φ m χ , while the exception is the s-channel resonance-enhanced region with m φ ∼ 2m χ . The latter case appears only for m φ (1) GeV to be consistent with collider constraints.

(b)
The lower limit on the WIMP mass is obtained by the combination of collider constraints and kinematic equilibrium condition. When it is lighter than (10) MeV, the mixing angle is constrained to be below 10 −3 due to Kaon experiments. The mixing angle is, however, required to be larger than 10 −3 to satisfy the kinematical equilibrium condition, as shown in the plot on the (m χ , | sin θ |)-plane in Fig. 8, which leads to the limit, m χ 10 MeV.
(c) On the other hand, the lower limit on m φ in the range of m χ 500 MeV comes from the ∆N eff constraint discussed in section 3.2.4, while the limit in the range of m χ 3 GeV is from the direct dark matter detection constraint, for the spin-independent scattering cross section between the dark matter and a nucleon is proportional to m −4 φ . The lower limit on m φ in the range of 500 MeV m χ 3 GeV is obtained in a complicated way: it is from the combination of the direct detection constraint and the kinematical equilibrium condition. When the dark matter mass is below a few GeV, the freeze-out temperature becomes lower than the QCD phase transition, so that a mixing angle is required not to be very suppressed in order to maintain the kinematical equilibrium. #22 However, such an unsuppressed mixing angle (as well as a small m φ ) leads to a large spin-independent scattering cross section between the dark matter and a nucleon and ruled out by the direct dark matter detection.
(d) In the right panel, the lower limit on the mixing angle is from the BBN constraint discussed in section 3.2.5, τ φ ≤ 1 s in the range of m φ ≥ 2m π and τ φ ≤ 10 5 s in the range of m φ < 2m π . The spike structure at m φ ∼ 1 GeV comes from the uncertainty on the decay width of the mediator discussed in section 2.4, while another one at m φ ∼ 200 MeV is from the quick change of the decay width due to the threshold of the φ → µ − µ + channel.
(e) On the other hand, the upper limit on the angle comes mainly from collider constraints. For m φ 500 MeV, the limit is set by Kaon experiments such as CHARM. For 500 MeV m φ 5 GeV, B meson experiments set the limit. The long island at m φ 1 GeV and 10 −3 | sin θ | 10 −1 is again due to the theoretical uncertainty of the φ-decay width. The reason why this island region is isolated is that the displaced vertex search of B → X s φ at the BaBar experiment ruled out the parameter region of 0.5 GeV m φ 1.5 GeV and 5 × 10 −4 | sin θ | 10 −3 . The step-like structure at m φ ∼ 2-4 GeV is because φ → g g, cc, τ − τ + opens and the branching fraction of the observable channel, e.g. φ → µ − µ + , is suppressed.
(f) When the mediator becomes heavier than 5 GeV, the limit on the mixing angle is set by the direct dark matter detection at underground experiments. This is because the dark matter is required to be heavier than the mediator in most of the parameter region, for the relic abundance condition is satisfied by the χχ → φφ annihilation process, and the direct dark matter detection put a sever limit when m χ 5 GeV. Only the exception is found in extended region at m φ ∼ several 10 GeV and | sin θ | ∼ 10 −2 , where the relic abundance condition is satisfied by the s-channel resonance-enhanced annihilation χχ → φ → ff .
(g) The SN1987A constraint [52] is over depicted on the (m φ , | sin θ |)-plane as a darker transparent color region at 5·10 −7 | sin θ | 3·10 −5 and m φ (100) MeV. It is seen that the SN1987A constraint is potentially important, for it restricts the parameter region that is not restricted by other constraints. The exclusion region by the SN1987A constraint can be understood as follows: The mediator interaction is too weak to affect the SN cooling when | sin θ | 5 · 10 −7 , while the mediator decays or trapped inside the SN when | sin θ | 3 · 10 −5 and does not contribute the cooling. When m φ T c 30 MeV with T c being the critical core temperature of SN1987A, the mediator cannot be thermally, thus efficiently, produced.

Future prospects
Results of our analysis for the future prospects of the minimal WIMP model projected on the (m χ , m φ )and (m χ , | sin θ |)-planes are shown in the left and right panels of Fig. 4, respectively, where we have assumed that dark matter signals as well as mediator signals are not detected at any near future experiments/observations discussed in the previous section. #22 The process φ f → φ f via the h − φ − φ coupling does not work (with ' f ' being a SM fermion), for only light SM fermions exist in the universe and their couplings to h are suppressed by small Yukawa couplings.
Our results projected on all planes of input parameters are again found in appendix C.2.
(h) In the left panel of Fig. 4, comparing with the one in Fig. 3, the lower limit on m φ becomes severer as m φ 20 MeV in the range of m χ 300 MeV due to the futureexpected constraint from the ∆N eff measurement. Together with the relic abundance condition, it also gives a lower limit on the WIMP mass. On the other hand, the lower limit on the mediator mass in the range of m χ 300 MeV becomes stronger than those in Fig. 3 because the significant upgrade will be expected at the direct dark matter detection in the near future. The void region, which is located at m χ ∼ 1 GeV and m φ ∼ a few hundred MeV, is because of the constraint from the SHiP experiment. In addition, the mixing angle in this region is required to be enough large from the kinematical equilibrium condition when m χ ∼ 1 GeV, while such a mixing angle will be ruled out if no signal is detected at the experiment.
(i) In the right panel of Fig. 4, the lower limit on the mixing angle | sin θ | in the range of m φ 10 GeV is not very much different from the previous one in Fig. 3. On the other hand, when m φ 10 GeV, the region of | sin θ | (10 −9 ), where it was survived in Fig. 3, is removed in Fig. 4. This is because the search for the Higgs invisible decay at the HL-LHC experiment rules out the region if no signal of the decay is detected there.

(j)
The excluded region at m φ 100 MeV and | sin θ | ∼ 10 −6 is from the direct dark matter detection at future underground experiments. This is because, as discussed in appendix C, the dark matter mass must be larger than 200 MeV to satisfy the kinematical equilibrium condition in this parameter region, and a small mediator mass leads to a large scattering cross section between the dark matter and a nucleon despite of a suppressed mixing angle. The sensitive future direct dark matter detection thus start proving this region.
(k) On the other hand, the upper limit on the mixing angle is also very much improved.
The SHiP experiment will change the landscape significantly at the region of 0.3 GeV m φ 4 GeV and 2 × 10 −6 | sin θ | 10 −4 , if no mediator signal is detected. Moreover, LHCb and Belle II experiments could also put severe constraints on | sin θ | in the range of 0.5 GeV m φ 5 GeV, while the direct dark matter detection at future underground experiments could put a severer constraint on the mixing angle | sin θ | in the range of m φ 5 GeV.
(l) Finally, let us comment on the resonant annihilation region, where the relation m φ ∼ 2m χ holds. As seen in the right panel of Fig. 4, this region could still survive in the near future at m φ ∼ 2m χ ∼ 10 GeV and | sin θ | ∼ 10 −2 , though it obviously shrunk compared to the region in Fig. 3. This is because that the sensitivity of the direct dark matter detection is significantly improved when the dark matter mass is greater than several GeV.

Implication of the results
Here, we consider implication of the results discussed in previous two sub-sections. In Fig. 5, some results projected on various observables are shown, where the region consistent with all present experimental results discussed so far is depicted using the same color code as those in Figs. 3 and 4. Moreover, the region which will survive (at 95 % C. L.) even if no  dark matter and mediator signals are detected in the near future is also over depicted as a red line. We discuss below implication of the results in each panel of Fig. 5 in some detail.
In the top left panel, the results are projected on the (m χ , σ SI )-plane, where σ SI is the spin-independent scattering cross section between the WIMP and a nucleon, as defined in section 3.2.3. As expected, the direct dark matter detection is effective to directly search for the WIMP as far as σ SI is large enough. #23 On the other hand, it can be also seen that the wide region with a very small value of σ SI remains survived even in future, especially when m χ is larger than (1) GeV. This originates in the survived parameter region with a suppressed mixing angle (sin θ ∼ 0), and thus other experiments are required to test it.
In the top right panel, we present the favored regions for cτ φ vs. Br(B → X s φ) with τ φ and Br(B → X s φ) being the lifetime of the mediator and the decay branching fraction of the B meson into a mediator and a hadronic system with a strangeness [135], respectively. It is found that the region with cτ φ (10 −3 -10 3 ) m will be well searched for thanks to the flavor experiments. Belle II and SHiP experiments will play significant roles for the search #23 The improvement of the lower limit on m χ in the near future is because of the provisional ∆N eff constraint. in regions of 10 −3 m cτ φ 1 m and 1 m cτ φ 10 3 m, respectively. On the other hand, the direct dark matter detection plays an important role to search for the region of cτ φ 10 −3 , for both the WIMP and the mediator are required to be heavy enough to have such a short mediator lifetime. Finally, we found that the direct dark matter detection and BBN constraints will play some roles to search for the region of cτ φ 10 5 m, though the wide parameter region survives even in the future due to the small mixing angle.
In the bottom left panel, we present the favored regions for cτ φ vs. Br(h → φφ), where Br(h → φφ) is the branching fraction of the exotic Higgs decay into two mediators. As discussed above, the direct dark matter detection will play a crucial role in the near future, in particular to search for the region of cτ φ 10 −5 m. On the other hand, other experiments will also play important roles to search for the region of cτ φ 10 −5 m. For instance, LHC experiment will search for the region by observing the exotic Higgs decay when its branching fraction is large enough, while cosmological observations (∆N eff measurement, etc.) and flavor experiments (Belle II and SHiP experiments, etc.) will search for the same region but when the branching fraction of the exotic Higgs decay is very suppressed.
It is worth emphasizing here that the constraint from the exotic Higgs decay plays a complemental role to others as shown in the bottom right panel. Magnitude of WIMP and mediator signals is proportional to the mixing angle squared in most of experiments and observations, while that of the exotic Higgs decay does not rely on the mixing angle. This is because the mediator is required to have an interaction to SM particles with enough magnitude to satisfy the kinematical equilibrium condition, and it is achieved by the φ−φ−h interaction (originating in the Φ 2 |H| 2 operator) when the mixing angle (originating in the Φ|H| 2 operator) is suppressed. As a result, the precise measurement of the nature of the Higgs boson will be mandatory in the future to test the minimal WIMP model, and it could be done by future lepton colliders such as ILC [150], CEPC [151] and FCC-ee [152].

Summary
We have studied a minimal (renormalizable) model with a light fermionic WIMP and a light scalar mediator whose decay width is computed with uncertainties from non-perturbative QCD effects properly taken into account. In order to investigate the present status and future prospects of the light WIMP and the light scalar mediator, we have performed a comprehensive likelihood analysis involving all robust constraints obtained so far. In addition, we also discuss those future sensitivities which will be obtained in the near future (if no WIMP signals are detected) from particle physics experiments as well as cosmological and astrophysical observations. We have carefully involved a kinematical equilibrium condition assuming that the (chemical) freeze-out of the light WIMP occurred when it was in kinematically equilibrium with the thermal bath composed of SM particles. We have paid particular attention to a possible case that the light WIMP can be in the kinematical equilibrium through existent mediator particles at the freeze-out epoch even if the WIMP does not have an interaction directly to SM particles with enough magnitude.
A very wide parameter region is still surviving at present in the region of 10 MeV m φ m χ , and it is found that many kinds of experiments and observations are required to test the WIMP in the near future. Direct dark matter detection experiments will play a crucial role in particular to search for the WIMP with the mass greater than (100) MeV, while flavor experiments will play a significant role to search for the mediator in the wide region of its mass. Moreover, precise cosmological observations such as the ∆N eff measurement are mandatory to search for a very light mediator whose mass is of the order of 10 MeV.
On the other hand, a wide parameter region will remain survived in the near future even if no WIMP/mediator signals are detected in the experiments. In principle, the experiments rely on the mixing angle (originating in the Φ|H| 2 operator) to detect WIMP/mediator signals, while all cosmological conditions (relic abundance and kinematical equilibrium conditions) can be satisfied even if the mixing angle is very suppressed. This is because the annihilation process within the dark sector, the χχ → φφ process, satisfies the relic abundance, while the kinematical equilibrium condition is satisfied through the other interaction between the mediator and SM particles, namely the φ −φ −h interaction (originating in the Φ 2 |H| 2 operator). As a result, a future experiment which is sensitive to this interaction is mandatory to test the remaining parameter region, and it can be achieved by measuring the nature of the Higgs boson precisely. In particular, future lepton colliders such as ILC [150], CEPC [151] and FCC-ee [152] will provide an ideal environment for this measurement.

A Preselection criteria
As mentioned in section 3.1, we apply preselection criteria on model parameters as c p = 0, |θ | ≤ π/6 and V (η, ξ) ≥ V (v Φ , v H ) in the range of |ξ| ≤ 1 TeV and |η| ≤ 1 TeV, where ξ and η are defined as Φ = η and H = (0, ξ/ 2) T , respectively. The region of the parameters survived after applying the criteria is shown in Fig. 6. The survived parameter region is further investigated by the subsequent likelihood analysis, as discussed in section 3.
Since the mediator mass is eventually constrained to be less than, or at least the same order of m χ to satisfy the relic abundance condition, as seen in section B.2, let us focus on the region of m φ (10) GeV and discuss how the preselection criteria work to restrict the model parameter space of the model. First, it can be seen in the result on the (m φ , µ 2 Φ )-plane that the parameter µ 2 Φ is constrained to be |µ 2 Φ | (100 GeV) 2 , because otherwise m φ would be much heavier than (10) GeV. Next, the quartic coupling of the field Φ is constrained to be λ Φ 0 as seen in the result on the (m φ , λ φ )-plane to make our vacuum stable, because the contribution from the quadratic term µ 2 Φ Φ 2 is small. Finally, the parameter of the cubic term is constrained to be |µ 3 | 100 GeV, as seen in the result on the (m φ , µ 3 )-plane, to suppress the negative contribution to the potential in the region of Φ = η 0.

B Kinematical equilibrium condition
In this appendix, we first write down all the reaction rates Γ χSM Γ χφ and Γ φSM in a concrete form. Then, we present the model parameter region survived after applying the kinematical equilibrium condition as well as the preselection criteria and the relic abundance condition.

B.1 The reaction rates
We first consider the reaction rate between the WIMP and SM particles Γ χSM . Considering the fact that the WIMP is always non-relativistic while SM particles are relativistic during the freeze-our process, the reaction rate Γ χSM is described by the following formula: where n SM i (T f ) is the number density of the SM particle 'i' at the freeze-out temperature.
Since relativistic SM particles dominantly contribute to the rate, its approximate form is given by Here, g i is the color and spin degrees of freedom of the particle 'i'. #24 The prefactor 2T f /m χ comes from the fact that about m χ /(2T f ) times collisions are required to maintain the kinematical equilibrium, because the typical momentum transfer in a single collision is estimated to be ∆q 2 /q 2 ∼ 2T f /m χ 1 in a stochastic process [153][154][155]. The thermal-averaged scattering cross section 〈σ χSM i v〉 T f is calculated according to the following formula [156]: where K n (z) is the modified Bessel function of the 2nd kind. When the freeze-out temperature is less than the QCD scale, namely T f ≤ Λ QCD 155 MeV [157], the scattering process between the WIMP and an electron or a muon by the exchange of φ in the t-channel dominantly contributes to the reaction rate. On the other hand, when T f ≥ Λ QCD , scattering processes with various SM particles contribute to the rate. Those are scattering processes of the WIMP with a muon, a tau lepton as well as strange, charm and bottom quarks.
We next consider the reaction rate between the WIMP and the mediator, Γ χφ . Since the mediator can be either relativistic or non-relativistic at the freeze-out temperature depending on its mass m φ , the rate is given by the formula which is similar to that in eq. (32): where n φ (T f ) is the number density of the mediator at the freeze-out temperature. The prefactor F (x 1 , x 2 ) is the extension of that in eq. (32), and it can be estimated based on the following discussion. First, when x 1 1 & x 2 1, or x 1 1 & x 2 1, the prefactor should be the same as the one in eq. (32), namely F (x 1 , x 2 ) = 2/x 1 or 2/x 2 . Next, when both arguments are small enough, the prefactor F (x 1 , x 2 ) should be (1), for the momentum #24 The exact formula of the number density, d 3 k g i /[exp(−E k /T ) ∓ 1], is used in our numerical analysis. transfer ∆q 2 /q 2 between relativistic particles is so large that a single collision is (almost) enough to maintain the kinematical equilibrium. Detailed calculation of the momentum transfer between relativistic particles expects that the prefactor is F (x 1 , x 2 ) 1/2. Finally, when the both arguments are (almost) the same and large enough, the prefactor should be (1) again, for the momentum transfer is expected to be efficient between two nonrelativistic particles whose masses are (almost) the same. Taking the above discussion into account, we use the following formula for estimating the prefactor F (x 1 , x 2 ): The scattering process between the WIMP and the mediator φ comes from diagrams exchanging the WIMP in the sand u-channels as well as a diagram exchanging the mediator in the t-channel. Its thermal-averaged scattering cross section, 〈σ χφ v〉, is calculated using the same formula as that in eq. (33) with m SM i being replaced by the mediator mass m φ .
We finally consider the reaction rate between the mediator φ and SM particles. Three different processes contribute to the rate: the (inverse) decay process Γ φ↔SMs , the scattering process Γ φSM↔φSM and the absorption (emmision) process Γ φSM↔SMs . The contribution from the (inverse) decay process is given by the thermal-average of the total decay width of the mediator Γ φ that has been discussed in section 2.4. Its explicit form is given as follows: where γ is the so-called the Lorentz gamma factor, and the distribution in the square brackets is the Maxwell-Juttner distribution. The contribution from the scattering process is given by a similar formula discussed in eq. (32). Since the mediator can be non-relativistic or relativistic at the freeze-out temperature T f , while relativistic SM particles dominantly contribute to the reaction rate, the explicit form of the contribution is given as follows: When T f ≤ Λ QCD , the scattering process between φ and an electron or a muon dominantly contributes to the reaction rate. On the other hand, when T f ≥ Λ QCD , scattering processes of φ with a muon, a tau lepton as well as strange, charm and bottom quarks contribute to the rate. Corresponding thermal-averaged scattering cross section are calculated using the same formula as that in eq. (33) with m χ being replaced by m φ . The contribution from the absorption (emmision) process is given by the same formula as that in eq. (37) again: where the scattering cross section σ φSM i in eq. (37) is replaced by the cross section of the absorption (emmision) process, σ φSM i . After the QCD phase transition, T f ≤ Λ QCD , the following three processes contribute to this reaction rate: where f is an electron or a muon, whilef is its anti-particle. On the other hand, when T f ≥ Λ QCD , in addition to the above processes, other three processes, φ g → qq, φq → gq and φq → gq also contribute to the rate with q (q) being a quark (anti-quark). Their corresponding thermal-averaged cross sections are calculated by the same formula as that in eq. (37). As a result, the reaction rate between the mediator φ and SM particles is given by the sum of all processes mentioned above: Γ φSM = Γ φ↔SMs + Γ φSM↔φSM + Γ φSM↔SMs .
As a demonstration, we show how each concrete process contributes to the reaction rates Γ χSM and Γ φSM in the left and the right panels of Fig. 7, respectively, at around the freeze-out temperature T f . Model parameters are set to be the same as those adopted in Fig. 2, namely (m χ , c s , m φ , sin θ , µ 3 ) are fixed to be (200 MeV, 0.022, 100 MeV, 10 −3 , 10 MeV).

B.2 Parameter region after the equilibrium condition applyied
Here, we present the model parameter region survived after applying the kinematical equilibrium condition as well as the preselection criteria and the relic abundance condition. As we mentioned in section 3.1.4, model parameters (m χ , c s , m φ , θ , µ 3 ) among seven independent parameters (m χ , c s , m φ , θ , µ 3 , µ 2 Φ , λ Φ ) are relevant to the following discussion, so that we show the result in Fig. 8 for the five parameters. In other words, parameters µ 2 Φ and λ Φ become nuisance parameters in the following phenomenological studies in this paper.
First, from the result on the (m χ , m φ )-plane, the mediator mass shall be at most around the WIMP mass to satisfy the relic abundance condition. It leads to the fact that the mediator is lighter than (10) GeV as long as we are discussing the light WIMP scenario. The upper limit on the mediator mass is also seen from the results on other planes spanned by m φ .      Third, the result on the planes spanned by the tri-linear coupling parameter µ 3 shows that the parameter is indeed restricted to be around (100) GeV, as addressed in appendix A.
Fourth, the coupling constant between the dark matter and the mediator |c s | is bounded from below due to the relic abundance condition, as can be seen from the result on the (m φ , |c s |)-plane. Here, one might worry about that the relic abundance of the WIMP becomes too small in the region with small m χ and large |c s |, however this region is satisfied by the case where m χ is slightly lighter than 2m φ but enough larger than m φ , namely the relic abundance condition is satisfied by the χχ → φ → ff process instead of χχ → φφ.
Fifth, the reason why the result on the (m φ , |c s |)-plane is similar to that on the (m χ , |c s |)plane is simply because the relation m φ ≤ (m χ ) holds in the whole parameter region.
Finally, on the (| sin θ |, |c s |)-plane, the reason why the region with small | sin θ | and small |c s | is excluded is as follows. When | sin θ | is small, m χ must be heavy enough to satisfy the kinematical equilibrium condition, as seen from the result on the (m χ , | sin θ |)-plane. On the other hand, m χ must be light enough when c s is small to satisfy the relic abundance condition, as seen in the result on the (m χ , |c s |)-plane. As a result, small | sin θ | and small |c s | are not simultaneously realized due to the contradiction between the two conditions.

C Supplemental figures C.1 Present status
Results of our analysis for the present status of the minimal WIMP model projected on all planes of input parameters (m χ , c s , m φ , θ , µ 3 ) are shown in Fig. 9. Since results on the (m χ , m φ )and (m φ , | sin θ |)-planes are already discussed in section 4.2, we will focus mainly on those on other planes, spanned by different combinations of the input parameters.
First, the result on the (m χ , | sin θ |)-plane can be understood from the discussion in sec-tion 4.2, because the relic abundance condition requires m χ m φ in most of parameter region. In fact, the lower limit on m χ , the lower limit on | sin θ | in the range of m χ 2 GeV and the upper limit on | sin θ | are understood in this manner. Only the exception is about the lower limit on | sin θ | in the range of m χ 2 GeV, where this region is excluded by the kinematical equilibrium condition, as already discussed in the previous appendix B.
Next, as seen from the result on the (m φ , µ 3 )-plane, the tri-linear coupling µ 3 is more constrained than the one in the corresponding panel of Fig. 8 when m φ 1 GeV. We have confirmed that this is due to constraints from meson decay experiments as well as the vacuum stability condition which was imposed as one of preselection criteria: The mixing angle | sin θ | is suppressed less than (10 −2 ) because of the collider constraints, so that µ 3 is required to be small enough when m φ 1 GeV in order to stabilize our vacuum. In a similar manner, results on the (m χ , µ 3 )and (| sin θ |, µ 3 )-planes can be understood as well.
Third, concerning the result on the (m φ , |c s |)-panel, the lower limit on m φ comes mainly from the ∆N eff constraint, though the upper-left corner with | sin θ | 5 × 10 −2 is further constrained by the direct dark matter detection. The lower bound on the coupling |c s | in the bulk region is from the relic abundance condition, where it is satisfied by the χχ → φφ annihilation. The result on the (m χ , |c s |)-plane is also understood in the same manner.
Fourth, the result on the (µ 3 , |c s |)-plane can be understood by those on the (m χ , |c s |)-and (m χ , µ 3 )-planes. When |c s | is smaller, a larger m χ is not allowed due to the relic abundance condition. On the other hand, a smaller m χ leads to a very restricted tri-linear coupling µ 3 .
Finally, the result on the (| sin θ |, |c s |)-plane can be understood as follows: The region of | sin θ | 10 −3 is not very much different from that in the corresponding panel of Fig. 8. On the other hand, the region of 10 −3 | sin θ | 10 −1 is excluded by meson decay experiments, because a smaller |c s | indicates a smaller m φ , as seen on the (m φ , |c s |)-plane. The region of | sin θ | 10 −1 is excluded by collider experiments as well as the direct dark matter detection.

C.2 Future prospects
Results of our analysis for the future prospects of the minimal WIMP model projected on all planes of input parameters are shown in Fig. 10, assuming that no dark matter and mediator signals are detected even in the near future. Since results on the (m χ , m φ )and (m φ , | sin θ |)planes are already discussed in section 4.3, we will focus on those on other planes spanned by different combinations of the input parameters, as in the previous subsection C. 1.
First, the result on the (m χ , | sin θ |)-plane can be understood from the discussion in section 4.3: the allowed region shrank compared to the corresponding panel in Fig. 9 because of the provisional future-update on ∆N eff , direct detection and collider constraints.
Next, the results on (m φ , µ 3 )-, (m χ , µ 3 )and (| sin θ |, µ 3 )-planes show that the constraint on the tri-linear coupling µ 3 becomes severer than those of Fig. 9 even at m φ m χ a few GeV, because future collider experiments will put a severer constraint on these mass region, if no signal is detected there. Moreover, µ 3 is highly restricted in the range of m φ a few ten GeV or | sin θ | 10 −3 , but the resonant annihilation region is severely constrained by the provisional direct dark matter detection in the near future, as discussed in section 4.3.
Third, allowed regions on (m χ , |c s |)-and (m φ , |c s |)-planes are, overall, shrunk compared to those in Fig. 9. Those lower limits on m φ and m χ come from the ∆N eff measurement, while the lower limit on |c s | is from the relic abundance condition. The small void region at 0.3 GeV m χ 1 GeV (0.3 GeV m φ 1 GeV) and |c s | 0.1 is due to the constraint from Moreover, the allowed region on the (µ 3 , |c s |)-plane remains similar as that of the present status in Fig. 9, except the resonant annihilation region at |c s | ∼ 10 −2 is shrinking.
Last but not least, on the (| sin θ |, |c s |)-plane, it is seen that the lower limit on the coupling constant |c s | is more or less the same as that in Fig. 9. It is again from the relic abundance condition as seen in the (m φ , |c s |)-and (m χ , |c s |)-planes. On the other hand, future meson decay experiments and direct dark matter detection make the parameter region of | sin θ | 10 −3 excluded. On the other hand, the shape of the contour at the region of | sin θ | 10 −9 is mainly due to the uncertainty of the mediator decay width as seen in the (m φ , | sin θ |)-plan. Since the dark matter mass is almost fixed to be around 10 GeV as seen in the (m χ , | sin θ |)plan, it requires a specific value of c s to satisfy the relic abundance condition.

D Relaxing Kinematic Equilibrium Condition
In the main text, we imposed the kinematic equilibrium condition at around the freeze-out temperature T D ∼ T f to figure out the conventional WIMP parameter region. However, it is also acceptable that the dark sector and the SM sector kinematically decoupled at higher temperature above the freeze-out, and then two sectors evolve independently [42].
Note that we adopt the condition T D = T f to figure out a very conventional WIMP parameter region in our setup. On the other hand, the condition can be relaxed by requiring that the WIMP is in the equilibrium at some temperature of the universe before the freezeout, because it still allows us to make a quantitative prediction on its abundance. In this section, we will discuss how the result of our analysis alters by relaxing the condition.
In order to understand how the kinematic equilibrium condition affects the parameter region, we relax the condition by requiring the decouple temperature T D above the freezeout temperature such as T D = 10T f or T D = 100T f . Then, we show how the thermal DM parameter region is expanded in the (m χ , | sin θ |)-and (| sin θ |, |c s |)-planes in Fig. 11.
The expansion of the parameter region can be understood as that the higher decoupling temperature T D increases the light degree of freedom from heavier SM particles, which have larger Yukawa couplings to help maintaining the kinematic equilibrium. From another point of view, because the freeze-out condition fixes the relation T f m χ /20, relaxing the decoupling temperature T D to higher temperature than T f allows a lighter DM mass region still maintaining the same degree of freedom and keeping the thermal equilibrium.
In the (m χ , | sin θ |)-plane in Fig. 11, for the T D = T f case, the vertical edge of m χ 2 GeV is from the pion threshold, which shifts to m χ 0.2 GeV and 0.02 GeV for the T D = 10T f and T D = 100T f cases, respectively. We can clearly see how the parameter region is extended to the one with a smaller value of m χ . The same behavior can be seen in the (| sin θ |, |c s |)-plane in Fig. 11, where the parameter region extends to a smaller value of |c s | due to the correlation between m χ and |c s | seen in the (m χ , |c s |)-plane in Fig. 8.
It is worth pointing out that, once the kinematic equilibrium condition is relaxed, the temperature of the dark sector (both dark matter and mediator) during the freeze-out could be different from the temperature of the SM thermal bath. Since the result in Fig. 11 is obtained assuming both the temperatures are (almost) equal, above discussions are validated only in such a case. #26 On the other hand, when the temperatures are very different, the relic abundance condition as well as BBN and N eff constraints have to be altered and the survived parameter space will be changed accordingly. A comprehensive study of early decoupled scenarios are indeed interesting but beyond the scope of our current study.