A tale of two dark neighbors: WIMP n’ axion

We study the experimental constraints on a model of a two-component dark matter, consisting of the QCD axion, and a scalar particle, both contributing to the dark matter relic abundance of the universe. The global Peccei-Quinn symmetry of the theory can be spontaneously broken down to a residual Z2-symmetry, thereby identifying this scalar as a stable weakly interacting massive particle, i.e., a dark matter candidate, in addition to the axion. We perform a comprehensive study of the model using the latest data from dark matter direct and indirect detection experiments, as well as new physics searches at the Large Hadron Collider. We find that although the model is mostly constrained by the dark matter detection experiments, it is still viable around a small region of the parameter space where the scalar dark matter is half as heavy as the Standard Model Higgs. In this allowed region, the bounds from these experiments are evaded due to a cancellation mechanism in the dark matter-Higgs coupling. The collider search results, however, are shown to impose weak bounds on the model. ar X iv :1 81 0. 09 47 1v 1 [ he pph ] 2 2 O ct 2 01 8

mass inversely proportional to the U (1) PQ -breaking scale. In the early universe, axion can be produced non-relativistically through a coherent oscillation of the axion field due to the misalignment of the PQ vacuum. This is known as the misalignment mechanism of axion production [21,22]. The axion is not completely stable, however, it has very feeble couplings with SM particles, thereby ensuring a lifetime longer than the age of the universe [23]. This makes the axion a very good CDM candidate, although the same feeble couplings make direct detection of these axions challenging [24].
In this work, we study a two-component DM model consisting of a WIMP and the axion as the DM candidates. As a simple realization of this, one can consider the KSVZ (Kim-Shifman-Vainshtein-Zakharov) model [19,20] of axion with an additional scalar field charged under the U (1) PQ [25]. This additional scalar gets its stability from the residual Z 2 -symmetry of the broken U (1) PQ , and hence becomes a WIMP-like DM candidate [26]. Breaking of the U (1) PQ and the electroweak symmetry leads to a mixing between the Higgs and the radial part of the PQ scalar, which leads to interesting phenomenological consequences. The advantage of this model is that although the axions have very weak interactions with the SM, the coupling between this dark scalar and the SM Higgs doublet provides a portal to test this model in different DM detection experiments, both direct and indirect. The model can also give different signatures at collider experiments. For example, the KSVZ model predicts new colored, electroweak singlet quarks, which can be produced at colliders. Mixing with a scalar affects the properties of the Higgs boson, which can be directly used to constrain the mixing parameters. Furthermore, the dark scalar can also contribute to momentum imbalance in a collision event.
Hence, in the light of recent experiments, we explore the constraints on the WIMPaxion DM model, both from DM search experiments, as well as collider searches. Using the recent limits on DM direct detection from XENON1T×1yr experiment data [27], we find that the phenomenologically interesting mass range of m DM 100 GeV is ruled out in such models. However, the stringent bounds from XENON1T×1yr data can be evaded in a small region of the parameter space where the scalar dark matter is half as heavy as the Higgs. This is a direct outcome of the mixing of the Higgs with the scalar, which leads to a cancellation mechanism in the Higgs portal coupling, thereby reducing the DMnucleon scattering cross-section. As a result, while minimal scalar DM models are mostly ruled out by direct detection bounds [28], such WIMP-axion models can still survive with a reduced parameter space. Collider signals, on the other hand, are highly plagued by the backgrounds from the production of standard model particles, and hence the signals are not significant enough to be observed above the background [29][30][31].
The paper is organized as follows. Section 2 discusses the model, and the different parameters involved. Section 3 talks about the different experimental bounds, and how they constrain the parameters of the model. In section 4, we summarize the main results, and finally in section 5, we conclude.

The Model
We consider the KSVZ model of the axion, where electroweak singlet quarks Q L , Q R and a complex scalar ζ, both transforming under a global U (1) PQ symmetry, are added to the SM [19,20]. These quarks are vector-like, hence do not introduce any chiral anomaly [32,33]. We augment this model with a complex scalar χ = (χ 1 + iχ 2 )/ √ 2 which is a SM singlet but charged under the U (1) PQ symmetry [25]. The axion a is the Nambu-Goldstone mode of the scalar field ζ, which can couple to the vector-like quarks, as well as χ. As in the original KSVZ model, the axion can act as a CDM candidate [24]. The charges and quantum numbers of the new particles are listed in table 1.
The relevant part of the Lagrangian, governing the interactions of Q L,R , ζ, and χ with the SM, is given by Here H is the SM Higgs doublet and d R represents right-handed down-type quarks in the SM. After electroweak symmetry breaking via the Higgs vacuum expectation value Similarly, using the nonlinear representation, one can write ζ = e ia /Fa (F a + σ 0 ) / √ 2, where F a is the U (1) PQ symmetry breaking scale as well as the axion decay constant, and σ 0 is the radial excitation of the ζ field. Constraints from supernova cooling data disfavor values of F a smaller than 10 10 GeV [34].
After the breaking of both the symmetries, viz. electroweak and PQ symmetries, the interaction term between H and ζ fields leads to mixing between h 0 and σ 0 with the mass matrix As a result of the mixing, the scalars in the mass basis are related to those in the flavor basis as where the mixing angle, in the limit F a v H , is given by One obtains the masses of the physical states as Note that the mass m h of the mixed state h is no longer 2λ H v 2 H , as predicted by the SM. Since h is the physical state, we fix m h at 125 GeV and the Higgs vev v H at 246 GeV to match with the experimentally measured masses of the observed scalar [35,36] and W, Z bosons respectively [37]. The value of λ H is no longer the SM value, λ SM H 0.13, but is dependent on other parameters in this model and can be calculated using Eq. (2.5).
, from Eq. (2.5) it is evident that λ ζH has to be zero, i.e., the SM Higgs does not mix with ζ, as considered in [25]. Note that there is no underlying symmetry in the theory that allows us to set λ ζH to zero in the Lagrangian. More importantly, although the mixing is very small, the relation of the masses of the physical state with other model parameters plays a major role in imposing constraints on the model. Therefore, we do not neglect the mixing of h 0 with σ 0 in this study. The masses of χ 1 and χ 2 are given by Without loss of generality, we can take χ > 0 such that m χ 1 < m χ 2 , hence χ 1 can be the DM candidate, and we, henceforth, denote the mass of χ 1 as just m χ . Note that after the PQ-symmetry breaking, the Lagrangian in Eq. (2.1) has a residual Z 2 -symmetry which stabilizes χ 1 . Note that in Eq. (2.7), µ 2 χ is defined to be negative and hence cancels out the large contribution coming from F a . This type of fine-tuning is a general feature of these axion models [25]. Since the fine-tuning is required mainly in the dark sector, we do not explore it further and defer the details to a later work. Furthermore, one can also motivate a tiny value of χ from naturalness arguments. As χ → 0, one obtains an extra U (1) symmetry in the theory, apart from the U (1) PQ . This can allow χ to be naturally small.
The mass of the axion is obtained through non-perturbative QCD effects and is inversely proportional to F a , m a 0.6 meV × 10 10 GeV F a . (2.8) The couplings of the axion to SM particles are also suppressed by inverse power of F a , so the decay lifetime of the axion is very large. In fact, if we take the value of F a > 10 10 GeV, as allowed by the supernova cooling data [34], its lifetime becomes larger than the age of the universe. Thus, the axion also acts as a viable candidate for CDM in this model. Therefore, both χ 1 and the axion will contribute to the total DM relic density in the universe. Finally, the vector-like quarks obtain their mass m Q = f Q F a / √ 2, as ζ develops a vev. If this mass is ∼ O(TeV), they can be produced at the LHC. This is expected to give direct constraints on this model, however, in order to have a mass of ∼ O(TeV), the coupling f Q needs to be extremely tiny ∼ O(10 −6 ).
The new interactions introduce two portals connecting the SM and the dark sector through the Higgs (via the hχ 1 χ 1 ) and the down-type quark (via the χ 1QL d R ). Of the two, the hχ 1 χ 1 interaction is the more important one and will play a key role in our analysis. The hχ 1 χ 1 coupling is given by Though sin θ is small, the first term cannot be ignored due to the large scale F a . Using the approximation for sin θ in Eq. (2.4), we obtain Note that in the presence of nonzero λ ζH , the hχ 1 χ 1 coupling vanishes at 11) as opposed to at λ χH = 0 as in [25]. This shift will play a crucial role in the analyses to follow. Using Eq. (2.5), λ ζ can be written in terms of m h , λ ζH , and λ H . This gives a family of solutions, satisfying Eq. (2.11). In figure 1, we show four contours of λ H in the λ ζH − λ χH plane for a given value of λ ζχ = 0.1. Any point on these hyperbolae satisfies Eq. (2.11), leading to vanishing hχ 1 χ 1 coupling. The benchmark point chosen for further analysis, λ ζH = 0.1, λ χH = 0.14 and λ H = 0.2, is shown as a black circle on the plot. One can, in principle, probe other values of λ H in this parameter space, and we do not show them here for clarity. However, one should not take λ H < λ SM H 0.13 since it leads to negative values of λ ζ , thereby making the potential for ζ unstable.
Finally, note from Eq. (2.6) that the mass of σ is proportional to the U (1) P Q -breaking scale F a . So if λ ζ ∼ O(1), σ becomes very heavy and decouples from the low energy theory. Therefore, for all practical purposes, σ does not play any significant role in present experiments. However, it is possible to have the mass of σ at around TeV, but only within a highly fine-tuned region of the parameter space.
One may wonder as to how much fine-tuning might be necessary for this scenario. Without going into details, we provide a back-of-the-envelope estimate here. From Eq.  However, in order to keep the physical masses real, i.e., both the eigenvalues of the mass matrix positive, the off-diagonal terms have to be of the same order as the diagonal terms. This requires λ ζH to be further fine-tuned to values ∼ 10 −7 . However, such small values of λ ζ and λ ζH will raise the value of g hχ 1 χ 1 (see Eqs. (2.9) and (2.10)) to values 1, which makes the whole problem highly non-perturbative. Then, one would again need to choose λ ζχ unnaturally small to solve this issue. 1 Since the above scenario is fine-tuned, we do not pursue it here. Rather, we consider natural values of all couplings O(1). As a result, in this work, the heavy scalar σ decouples early on and does not enter our analysis.
2. Mixing between h 0 and σ 0 changes the couplings of the observed 125 GeV scalar from that of the SM Higgs. This leads to changes in the properties of the observed scalar measured in the collider experiments from that of SM Higgs. This will also constrain the parameters of the model. 3. Since the masses of the DM and the vector-like quarks are lighter or near TeV range, they can potentially be produced at the LHC. Non-observation of such particles will limit the model parameter space.
The rest of the section discusses these types of experimental constraints in details.

Dark Matter Relic Abundance
After the U (1) PQ -symmetry breaking, the axion a, being a Nambu-Goldstone, enjoys a continuous shift symmetry. This symmetry is broken explicitly as a result of the chiral symmetry breaking in the QCD sector, and a temperature-dependent potential for the axion is generated from non-perturbative QCD effects [40]. But the axion field does not start rolling in the potential and remains frozen at its initial value until its mass becomes larger than the Hubble expansion rate is the scale factor of the universe. After the epoch when m a (t) H(t), the field starts oscillating coherently and the axion particles are produced with non-relativistic speed. They contribute towards the CDM abundance today and their density is approximately given by [24,41], Here θ a is the initial misalignment angle of the axion field relative to the minimum of the axion potential. For simplicity, we shall assume θ a ∼ 1 in the rest of the paper [42]. In order that the axions do not overproduce DM in the universe, the PQ breaking scale F a has to be less than 10 12 GeV. In this work, we will focus on 10 10 GeV ≤ F a ≤ 10 12 GeV. As already noted, χ 1 gains stability from the residual Z 2 -symmetry and is a DM candidate. In the early universe, χ 1,2 are in chemical equilibrium with the thermal bath of the SM particles. As the temperature of the universe decreases below ∼ m χ /20, their rate of interaction drops below the expansion rate and χ 1,2 cease being in equilibrium with the SM particles. The heavier component χ 2 , however, does not remain stable since it decays to χ 1 , which then forms the relic abundance Ω χ h 2 . The relic abundance is formed after the freeze-out of χ 1 χ 1 annihilations. The annihilation can be mediated by h as well as σ. However, the h-mediated process dominates, since m σ m h . The relic abundance, being governed by χ 1 χ 1 → SM SM, depends directly on m χ .
We show the dependence of the χ 1 relic density as a function of its mass m χ in figure 2. We used micrOMEGAs5.0 [43] to numerically compute Ω χ h 2 . The behavior for very small and large m χ can be understood as follows. For very small values of m χ ( few GeV), χ 1 can annihilate only into light quarks and the cross section is suppressed by the small Yukawa couplings resulting in overabundance of χ 1 . For m χ m t , the annihilation cross section is 1/m 2 χ suppressed. Since the relic abundance is inversely proportional to the GeV is due to opening up of the χ 1 χ 1 → hh channel. The shaded region above the Ω c h 2 = 0.12 line is ruled out by the Planck experiment [1]. We allow the under-abundance regions as the axion may account for the rest of the relic abundance. Other parameters chosen for this plot are as follows: F a = 10 10 GeV, annihilation cross section, we expect the region around m χ ≈ 100 GeV to give the correct ballpark value of the desired relic abundance.
The sharp dip at m χ m h /2 62.5 GeV is due to the s-channel resonance from the h propagator. As m χ increases further from 62.5 GeV, the cross section falls leading to sharp increase in the relic. When the χ 1 is heavier than h, the new annihilation channel χ 1 χ 1 → hh opens up and dominates over all other channels. As a result, the relic abundance decreases, leading to the second dip. As χ 1 becomes more massive, the relic increases again because of the decrease in annihilation cross section with the characteristic 1/m 2 χ suppression. Note that we do not consider m χ > M Q , since the colored Q L,R can become the lightest dark sector particle.
In our analysis, we take the Planck (TT, TE, EE, lowP) measurement of the CDM energy density Ω c h 2 = 0.12 ± 0.0012 represented by the horizontal line labeled as Ω c h 2 in figure 2 [1]. The over-abundance region, shown as a gray shade, is disallowed. However, the under-abundance region is allowed since the axion abundance Ω a h 2 can account for the rest of the relic. Therefore, the observed relic abundance Ω c h 2 We note that Ω χ is virtually independent of F a due to v H /F a suppression in the couplings and mixing angle. Hence, F a is fixed by Eq. (3.2) via the Ω a h 2 term.

Direct Detection of Dark Matter Particles
The DM direct detection (DD) experiments look for scattering between the DM particle and nuclei of the detector material. Any interaction between the DM and the SM quarks/gluons in a given model leads to a possible signal in the direct detection experiments. Nonobservation of such a scattering signal in such experiments constrains the parameters of the model. In the present case, the dominant channel of interaction arises again through the hχ 1 χ 1 coupling, since h mediates the DM and SM quark scatterings. A typical behavior of the scattering cross section as a function of λ χH is shown on the left panel of figure 3. The cross section σ χN is constant for very small λ χH because the coupling becomes independent of λ χH . For very large λ χH , the cross section increases as ∼ λ 2 χH , as expected. In between, a dip occurs because of the cancellation of two terms appearing in the vertex factor of hχ 1 χ 1 coupling (see Eq. (2.10)). Presently, the most stringent bound on this cross section is given by the XENON1T×1yr experiment [27]. It is most sensitive to the DM mass in the range 10 GeV − 1 TeV and the strongest upper bound quoted is σ χN 10 −46 cm 2 . We will show later that due to the stringent constraint, the only experimentally allowed region of DM mass turns out to be around m χ 62.5 GeV.
The right panel of figure 3 shows the χ 1 -nucleon scattering cross section σ χN as a function of m χ for two different values of λ χH . Note that in this model, χ 1 forms only a fraction f χ (≡ Ω χ /Ω c ) of the present dark matter abundance. Therefore the XENON1T bound is to be accordingly divided by f χ before applying to this model.
All the above bounds apply for χ 1 as the DM candidate. However, direct detection experiments for axion need to follow a different search strategy because of its ultra-low  mass. There have been a few experimental efforts to look for axionic dark matter. For example, the ADMX experiment [44] uses RF cavity to look for its interaction with the electromagnetic field. In the KSVZ model, this interaction strength is given by [19,20] where α is the fine structure constant. Presently, ADMX rules out a narrow region of the parameter space above g aγ 10 −15 GeV −1 (F a 10 12 GeV) around m a 2 µeV. For higher mass axion, the bound is even weaker. Another proposed experiment is CASPEr-Electric which will probe F a 10 12 GeV for lighter axions [45]. Moreover, we should remember that these bounds assume that 100% CDM abundance is given by axion which is not be true in our model. These bounds are weaker than the upper limit on F a from the dark matter relic abundance, even after adjusting for correct factor to cancel out the assumption, hence does not need a special attention.

Dark Matter Annihilation Signal
Various astrophysical observations hint that the present day universe consists of galaxies sitting inside halo-like structures formed by gravitational clustering of DM particles [47]. At the center of these halos, the DM density is high enough to scatter with each other and annihilate into SM particles. These final state particles would further decay and give rise to gamma-ray signals from various astrophysical objects, such as dwarf galaxies, the Milky Way center etc. We focus on bounds arising from gamma-ray signals due to such annihilations of DM particles.
We pay more attention to the DM mass around m χ m h /2 = 62.5 GeV which is still allowed by the direct detection experiment data. The total annihilation is dominated by the bb-channel (∼ 90%), which is shown in figure 4. Note that here also the annihilation cross section is enhanced due to the s-channel resonance from the SM Higgs propagator. Hence the largest annihilation signal is predicted at this mass. The dependence on λ χH comes through the g hχ 1 χ 1 coupling. The Fermi-LAT constraint becomes ineffective for the value of λ χH for which this coupling vanishes (the red curve in figure 4) as is discussed in section 2.
There have been many experiments which have looked for DM annihilation signals from various astrophysical objects [46,[48][49][50]. At present, the most stringent upper bounds on the thermally averaged DM annihilation cross section σv is given by the DES-Fermi-LAT joint gamma-ray search from the satellite galaxies of the Milky Way [46]. It is derived from 6 years observation of 45 such objects by the LAT. They have relatively less amount of visible baryonic matter and the DM population is expected to dominate their matter density. In figure 4, we show this upper bound on the annihilation cross section as the gray shaded region. This does not rule out most part of our parameter space, except a region of m χ around Higgs mass. In passing we also note that the DM mass needed for the resonantly enhanced annihilation signal in the bb-channel matches the result of the galactic center excess analysis done in ref. [51] within 1σ C.L. (also see [52]).

New Physics Searches at the LHC
In this subsection, we will focus on various signatures of the model at the LHC. The model has an extended scalar sector: apart from the SM Higgs boson h 0 , there exists a scalar DM candidate χ 1 and its heavier counterpart χ 2 , and another scalar field σ 0 , which is the radial component of ζ. As discussed earlier, h 0 and σ 0 mixes with each other giving rise to physical states h and σ. The mixing between σ 0 and h 0 changes the properties of h from that of the SM Higgs via its coupling to SM particles as well as to the new states present in this model. Since various properties of the observed scalar particle at the LHC resemble that of the SM Higgs boson, we expect some constraints on the parameter space of the model from the measurement of the properties of the observed 125 GeV scalar.
One of the measurements that provides relevant information about the properties of the observed 125 GeV scalar is its signal strength. If the scalar decays to X ∈ { ± , q, g, Z, W } and its conjugate,X, its signal strength is defined as  partial decay width of the new decay modes of h is Γ new , the signal strength of h decaying to any SM particle pairs XX can be written as where Γ tot SM is the total decay width of SM Higgs boson. In table 2, we tabulate the recent measurements of signal strength of the observed scalar h by both ATLAS and CMS collaboration at 13 TeV with ∼36 pb −1 integrated luminosity in different decay channels of h. The superscripts in the µ XX represent the production mode of the scalar h. For our analysis, we constrain the parameter space by imposing the value to be at 95% C.L. of the measured values, i.e., with ±2σ around the measured central value. Since, in the model, µ XX is always below unity, it is the lower bound at 95% C.L. which will actually put constraints on the parameters.
In the left panel of figure 5, we show the variation of the signal strength of h in W W * channel as a function of λ χH for two different masses of χ 1 . As expected from Eq. (3.5), the variation is a Lorentzian, with a narrow width governed by Γ tot SM and m χ . Since the coupling for h to χ 1 χ 1 , as given in Eq. (2.10), vanishes at λ χH = 2λ ζχ λ H − λ SM H /λ ζH (≈ 0.14 for the chosen benchmark point), the decay mode for the h vanishes at that point, and hence the µ XX becomes 1 around that point. The gray (green) shaded region shows the area disallowed at 95% C.L. by the measurements by CMS (ATLAS) as indicated in the plot, and the allowed region is shown in white. Although the measurements for different decay channels of h are listed in table 2 for completeness, we only plotted µ (ggF) W W * , which gives the strongest bounds from the signal strength measurement.
We also study the bounds from the invisible decay of h which arises from the decay  channel h → χ 1 χ 1 for m χ < m h /2 in this model. The BR of the decay can be written as The dependence of BR(h → χ 1 χ 1 ) with the parameter λ χH is plotted in the right panel of figure 5 for two different masses of χ 1 . As in the case with the signal strength, the BR(h → χ 1 χ 1 ) vanishes at the point where the coupling of h to χ 1 χ 1 , given by g hχ 1 χ 1 (see Eq. (2.10)), goes to zero. This feature is evident from the plot in the right panel of figure 5. Away from this point, the BR increases in both sides, tending to unity for high value of g hχ 1 χ 1 , which indicates that Γ new is the dominant decay mode, and all other modes are suppressed. Non-observation of this decay mode of the observed 125 GeV scalar at the LHC, therefore, places upper limit on the invisible decays of h. These upper limits are tabulated in table 3. In the right panel of figure 5, the gray (green) shaded region is the area disallowed at 95% C.L. by CMS [64] (ATLAS [63]) measurements on the invisible decays of 125 GeV scalar. It is therefore clear that only a small range of λ χH , for which the BR curves fall within the white region, is allowed by current measurements.
At this point, it is worth mentioning that the trilinear coupling of h is also modified due to the mixing with σ 0 , which will change the di-Higgs production rate. Measurements for the trilinear coupling of h as well as di-Higgs production have been carried out by both ATLAS [65] and CMS [66] in the di-Higgs channel. However, the upper bounds are well above the SM prediction due to lack of signal in the di-Higgs channel. Hence, much of the parameter space, especially the region of interest, of the model is not constrained by the measurement of trilinear coupling of h.
The model also predicts new particles at around GeV-TeV range, which can potentially be observed in a TeV collider. One such particle is the DM candidate, χ 1 , which is weakly interacting and does not decay within the detector. If it is produced in the collider, it will not be detected and will contribute to the missing momentum in an event. The other particles, within the observable range of TeV collider, are the vector-like quarks Q L and Q R . Since these quarks are colored, they can be produced in a hadron collider and subsequently decay to a down-type quark and a χ 1 . Presence of χ 1 will again contribute to the missing energy in the detector. Lack of agreement of such signals with those predicted at the TeV colliders will also put bounds on the parameter space of the model in consideration. Now, we turn to the discussion of direct production of the new particles at the LHC. The new particles, being charged under a PQ symmetry, should be produced in pairs. There are three different pairs of new particles that can be directly produced: QQ, Qχ 1 , andQχ 1 . Hence, these processes will contribute to the following final states: dijet (2j)+MET in case of QQ production, and monojet (j)+MET final state in case of Qχ 1 andQχ 1 production, where MET stands for missing transverse energy. In the rest of this section, we will discuss the constraints on the parameter space in view of the observation of the above-mentioned final states at collider.
Since the Qs are colored, the cross section for the production of QQ will be similar  to that of the SM quarks and will be suppressed for higher masses. Figure 6 shows the variation of total production cross section for QQ (in red) and for Qχ 1 andQχ 1 (in blue) in 2j+MET and in j+MET channels respectively at the LHC at 13 TeV. The production cross section of QQ in 2j+MET channel have negligible dependence on f d,s,b since the dominant parton-level process for the production is gg → QQ, which is independent of f d,s,b . Hence, the two red curves, solid for f d,s,b = 0.1 and dashed for f d,s,b = 1 coincide with each other. However, the cross section for Qχ 1 andQχ 1 in j+MET channels scales as f 2 d,s,b since the parton level process involved in the production is gq, gq → Qχ 1 ,Q χ 1 , whose amplitude is proportional to f d,s,b . Note that the only possible decay mode of Q is to a down-type quark and a χ 1 .
To estimate the signature of our model in collider experiments, events have been gen-erated at partonic level using MadGraph5 [67] with NNPDF2.3LO parton distribution function [68] using the UFO files generated by FeynRules [38,39] at center-of-mass energy of 13 TeV; partons in the final state have been showered and hadronized using the parton shower in PYTHIA 8.210 [69] with 4C tune [70]. Stable particles have been clustered into anti-kT [71] jets of size 0.4 (used by both ATLAS and CMS) using FastJet [72] software package; only the jets with P T more than 30 GeV have been considered for further analysis. In figure 7, we present some important and representative differential distribution of some observables as are considered by experimental collaborations to search for signals. | p T j |, which is the scalar sum of p T of all the jets. The major sources of the SM backgrounds for jets+MET are from the production of Z decaying to νν and W decaying to τ ν τ in events with jets. Also QCD events are potential sources to contribute to the same final state. The distribution for these three backgrounds are plotted in four panels of figure 7. SM background samples have been generated with at leading order (LO) using MadGraph5 [67] with NNPDF2.3LO parton distribution function [68] at center-of-mass energy of 13 TeV and PYTHIA 8.210 [69], with the same 4C tune [70] as used for generation of the signal sample, has been used for the simulation of fragmentation, parton shower, hadronisation and underlying event. The distribution for QCD, W +jets, and Z+jets backgrounds are plotted with gray, purple, and green respectively with the same color convention in all the four panels. From the figure, it is quite clear that the bumps for signals will not be significant enough to be observed above the expected fluctuation of the background.
Following the distribution in the experimental references [31,[73][74][75][76][77][78][79][80], we carried out our analysis with the same distribution. As discussed earlier, the direct production of new particles will contribute to 2j+MET and j+MET signals. There are few dedicated search in these channels to search for dark matter signals [31,[73][74][75]. Few other models, especially SUSY in R-parity conserving scenario, also lead to these kinds of signals. These searches have also been done by both CMS [76,77] and ATLAS [78][79][80]. Though the results are given in terms of SUSY parameters or effective theory parameters, one can recast the result for a given model and check for its consistency. But these searches do not yield any further constraint in the parameter space in the model. A dedicated search for this model may give a stronger constraint, but the analysis of such search is beyond the scope of this work.

Results
Our main results are summarized in figure 8. The relevant bounds coming from the different experiments are imposed on the region satisfying the DM relic density in the λ χH − m χ plane. The gray shaded region is ruled out by the relic constraints. We allow for both χ 1 as well as the axion to contribute to the DM relic density. Hence the white region, corresponding to the 2σ bound Ω c h 2 < 0.12, represents the allowed parameter space,  Figure 8. Allowed regions in the parameter space for the two-component axion-WIMP DM model. The gray shaded region shows the area ruled out by DM relic abundance constraint corresponding to the 2σ bound Ω c h 2 < 0.12 [1]. The black hatched lines show the regions of parameter space ruled out by the DM direct detection bounds from XENON1T×1 yr experiment [27]. The hatched region within the red curve is ruled out by the DM annihilation data from DES-Fermi-LAT experiment [46]. The blue shaded region show the bounds imposed due to the invisible decay modes of the Higgs, which is roughly 25% of its branching ratio [53,54]. The bound coming from the signal strength of the Higgs is shown in orange [63,64]. The white, unshaded region represents the allowed parameter space in this model. satisfying the relic density. As explained before, near m χ ≈ m h /2, the DM annihilation cross section is enhanced from the Higgs resonance, thereby decreasing the relic density of DM. This explains why the allowed region from relic is centered around m χ = m h /2. Furthermore, there is a particular set of parameters for which hχ 1 χ 1 coupling vanishes, leading to a rise in the relic density. This accounts for the peak-like structure in figure 8, which occurs at λ χH ∼ 0.14 for our choice of parameters. The black hatched lines show the regions of parameter space ruled out by the direct detection bounds from XENON1T×1 yr experiment. The hatched region within the red curve is ruled out by DES-Fermi-LAT joint gamma-ray search data from the Milky Way satellite galaxies. As is clearly seen, most of the allowed regions are ruled out, leaving behind a tiny window around in the m χ − λ χH plane. Clearly, this window is centered around m χ ≈ m h /2 and the value of λ χH for which the hχ 1 χ 1 coupling vanishes.
The blue shaded region shows the bounds imposed due to the invisible decay modes of the Higgs, which is roughly 25% of its branching ratio. More stringent bounds are imposed from the signal strength of the Higgs, which is shown in orange. These also help to rule out extra regions of the parameter space for larger as well as smaller values of λ χH . We have also checked that the LHC bounds from production of QQ are relatively weak, hence they do not impose any extra constraint on the model. Thus, from the above figure, one concludes that only a small fraction of the model can still be accommodated from existing experimental bounds. This region, however, enjoys the advantage of an accidental cancellation of the couplings near m h /2, thereby making it extremely difficult to rule out experimentally. This tiny window provides a breathing space for the model to survive.

Summary and Discussions
In this paper, we have performed a comprehensive study of a two-component dark matter model, consisting of the QCD axion, and an electromagnetic charge neutral scalar particle, both contributing to the relic density. The theory is symmetric under a global Peccei-Quinn symmetry, which can be spontaneously broken down to a residual Z 2 symmetry. For concreteness, we have considered a specific model: the KSVZ model of the axion, augmented with an additional complex scalar. After spontaneous breaking of the PQ symmetry, the residual Z 2 symmetry allows the lightest component of the complex scalar to be a DM candidate, apart from the axion. We have tested the model in the light of recent data from DM direct and indirect search experiments. Furthermore, we have also studied the different collider signatures of this model.
Although the observational and experimental constraints are found to be very restrictive, a synergy of the enhancement of DM annihilation from the Higgs resonance, and the vanishing of the coupling between the Higgs and the dark matter leave room for future experimental investigation of this model. A large portion of the parameter space predicts overabundance of χ 1 in the universe and hence is not viable. In the remaining underabundant region of χ 1 , the axion can form the dominant part of the CDM. The viability of the axion being the CDM is being tested in several ongoing experiments. The latest dark matter direct and indirect detection experiments results further constrain this model. Moreover, these results are expected to improve the bounds by few orders of magnitudes over the next few years which will subject this model to even tighter constraints. Although the bounds from the measurements of the properties of the Higgs at collider experiments are relatively weak, they still help to rule out an additional part of the parameter space. Future measurements of vector-like quarks at high luminosity and high energy operating modes of the LHC can shed further light on the viability of this model. Nevertheless, it is possible to add new particles to this simplistic model, e.g., an additional scalar, to enrich its phenomenology and evade some of the experimental bounds. This leaves room for future scopes of model-building and investigation of observable signatures in high energy experiments. In this work, we have calculated the prediction from our model with some natural choice for the couplings to compare with experimental data; but other values of the couplings can be explored to test the validity of the model on the basis of available experimental results. In conclusion, the two-component dark matter model, consisting of the WIMP and the axion, continues to survive, in spite of being tightly constrained.