Collider and dark matter searches in the inert doublet model from Peccei-Quinn symmetry

Weakly Interacting Massive Particles (WIMPs) and axions are arguably the most compelling dark matter candidates in the literature. Could they coexist as dark matter particles? More importantly, can they be incorporated in a well motivated framework in agreement with experimental data? In this work, we show that this two component dark matter can be realized in the Inert Doublet Model in an elegant and natural manner by virtue of the spontaneous breaking of a Peccei-Quinn U(1)P Q symmetry into a residual ℤ2\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ {\mathbb{Z}}_2 $$\end{document} symmetry. The WIMP stability is guaranteed by the ℤ2\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ {\mathbb{Z}}_2 $$\end{document} symmetry and a new dark matter component, the axion, arises. There are two interesting outcomes: (i) vector-like quarks needed to implement the Peccei-Quinn symmetry in the model may act as a portal between the dark sector and the SM fields with a supersymmetry-type phenomenology at colliders; (ii) two-component Inert Doublet Model re-opens the phenomenologically interesting 100-500 GeV mass region. We show that the model can successfully realize a two component dark matter framework and at the same time avoid low and high energy physics constraints such as monojet and dijet plus missing energy, as well as indirect and direct dark matter detection bounds.


Introduction
The nature of dark matter (DM) is one of the most puzzling mysteries of astroparticle physics. Dark matter which accounts for 27% of the total energy density of the Universe is often interpreted in terms of WIMPs. More precisely, only one species makes up the entire DM of the universe. Although, this minimal scenario having just one particle making up the entire DM sounds appealing, there is no observational evidence supporting this idea. In fact, both matter and radiation components of the Universe energy budget is comprised of more than one particle. Thus it is rather possible that DM is constituted by more than one particle.
Albeit, having as a second WIMP merely to account for some fraction of the relic density is not so compelling, unless two solid signals are observed consistent with the two WIMPs scenario. Hence, an ideal scenario would occur if the second DM component solves a major problem in particle physics such as the Strong CP problem. This is precisely our case understudy, where the second component is the axion field. As we shall see, the addition of axions to the WIMP paradigm has two motivating outcomes: the WIMP is JHEP10(2016)015 stability is naturally addressed; the strong CP problem in the Standard Model (SM) is solved. To put this framework into perspective we need to briefly review how this come into play.
WIMPs that are the most popular candidate for DM suggest a connection between DM physics and the weak scale. The stability of the WIMP is usually assumed to be due to the presence of ad hoc discrete global symmetry, such as a Z 2 symmetry, which prevents its decay. Another candidate is the axion [1,2], which is the pseudo Nambu-Goldstone of the breakdown of the U(1) P Q Peccei-Quinn (PQ) symmetry proposed to solve the strong CP problem [3] (see refs. [4][5][6] for a review). Under the assumption that the U(1) P Q symmetry is broken at an energy scale much higher than the electroweak scale, the axion can be an ultralight particle with faint interactions with all other particles [7][8][9][10], and allowed to have a lifetime larger than the age of the Universe. The axion contribution to the total DM energy density in the Universe also depends on the energy scale in which the U(1) P Q symmetry is broken [11].
Thus, the scenario in which both WIMP and axion make up the DM of the Universe is a natural and compelling framework. With that in mind we add a new and well motivated ingredient, the axion, in one of the simplest extensions of the Standard Model with a WIMP: the Inert Doublet Model (IDM), which contains an additional SU(2) L Higgs doublet with the lightest component stabilized by an ad hoc Z 2 symmetry [12][13][14].
In other words, we propose the axion as the DM companion to the IDM component H 0 . To this end we have developed a model based on the observation made in [15], where a U(1) P Q symmetry broken spontaneously into a Z 2 symmetry was advocated to stabilize the WIMP. 1 We tacitly assume that the U(1) P Q symmetry is protected against gravitational effects -which generate Planck-scale-suppressed symmetry breaking operators -by some sort of discrete symmetry (as in e.g. [23,24]) to avoid destabilization of the solution to the strong CP problem, and also of the WIMP [25]. The use of this global symmetry to stabilize the WIMPs is safe from gravitational effects which might violate the U(1) P Q [25], since only Planck suppressed operators of dimension six are present. To complete this two component DM system, at least a scalar singlet field hosting the axion a and a vector-like quark D are necessary in addition to the inert Higgs doublet whose lightest neutral component is the heavy DM [26][27][28]. The vector-like quark allows a simple implementation of the U(1) P Q symmetry, as in the KSVZ axion model [7,8], and acts as a portal connecting the SM and the dark sector. As a consequence of the residual Z 2 symmetry, the heavy vector quarks decay only to new heavy scalars and SM quarks, mimicking the phenomenology of R-parity conserving supersymmetry (SUSY) at colliders, including the classic SUSY signal of jets plus large missing energy.
As there is currently many experimental constraints on supersymmetry from the LHC searches, we performed, prior to the study of the multi-component DM scenario of our model, an investigation of the limits from the searches of jets plus missing energy and monojets at the LHC. After that, we focused on the main goal of the paper, which is JHEP10(2016)015 the study of our axion-WIMP DM scenario, pointing out the differences in relation to the typical IDM. The main finding is that, in contrast with the one-component DM in the IDM, the phenomenologically important mass interval 100 GeV ≤ M H 0 ≤ 500 GeV is re-opened, with the axion filling the role of the remaining DM.
Thus in summary, our work extends previous works in several key aspects; • We add the axion to the WIMP paradigm in the scope of the Inert Doublet Model, naturally addressing the WIMP stability (section 2).
• We investigate the collider phenomenology of the model (section 3).
• We show that one can successfully have a two component dark matter model with a wealth of experimental constraints (section 4).

The model
The model consists on a KSVZ type axion model [7,8] with an inert doublet H D , whose the lightest neutral component is stabilized by a residual Z D 2 symmetry that remains unbroken from the original PQ symmetry. Therefore, we will have two candidates for DM: the ultralight axion and the WIMP-like lightest component of H D .
The simplest way to implement the breaking U(1) PQ → Z D 2 is to break the PQ symmetry by a vev of a singlet scalar S of PQ(S) = 2 while all other fields carry integer PQ charges. The fields carrying even or zero PQ charge will be even under the remaining Z D 2 whereas those carrying odd PQ charge will be odd under Z D 2 , and thus belong to the dark sector. The conservation of Z D 2 requires that scalars with odd PQ charge be inert. As usual, the responsible for PQ symmetry breaking will host the axion in its phase as where a(x) is the axion field, and f a the axion decay constant that corresponds to the vev of S in our case (a KSVZ type axion model [7,8]). Nonperturbative QCD effects lead to a potential, which generates a mass to the axion as m a ≈ 6 meV × (10 9 GeV/f a ) .

(2.2)
In this framework the axion couplings with matter and gauge bosons are suppressed by f a which, being much higher than the electroweak scale, makes the axion an ultralight particle with feeble couplings to all other particles. In fact, f a is constrained from astrophysical objects which would have their dynamics affected if axions interact too much with photons. For example, supernova SN1987A data constrains f a to be greater than 10 9 GeV [29,30]. Still, an upper limit on the decay constant is obtained from the requirement that the axion relic density should not exceed the DM density, which gives f a ≤ 10 12 GeV [31][32][33][34].
In addition to the SM fermions we assume that there is at least one heavy quark field D ∼ (1, −1/3), where the numbers inside the parenthesis represent the transformation properties under the electroweak gauge group factors SU(2) L and U(1) Y ; the case of charge JHEP10(2016)015 Table 1. Quantum numbers of the fields beyond the SM.
2/3 exotic quark can be treated analogously. Such a quark field is formed by left-and righthanded fields D L,R , having the following interaction with S With the vev of S a mass M D = y f a / √ 2 for the D quark is generated through the interaction in eq. (2.3). We tune the Yukawa coupling y ≤ 10 −6 in eq. (2.3), as f a ≥ 10 9 GeV, so that M D lies at the TeV scale. In the appendix A.2 it is shown how to ameliorate such a tuning by extending the model.
Besides the fields necessary to solve the strong CP problem, we augment the SM with an inert Higgs doublet H D ∼ (1, 1/2), with PQ(H D ) = −1, in addition to the usual Higgs doublet H ∼ (1, 1/2). In the limit of PQ symmetry conservation, the Higgs potential is effectively 2 Exact PQ symmetry at the electroweak scale would imply degenerate CP odd and CP even scalars of the inert doublet, a feature that is problematic if the inert doublet accounts for all or most of the DM: direct detection searches for inelastic DM requires a mass splitting larger than 100 keV [28,35]. As PQ symmetry is broken at the scale f a we expect the additional PQ-violating but Z D 2 conserving term to be generated: The mass splitting is thus controlled by λ 5 , which can be taken real. A simple completion that generates the term (2.5) is shown in appendix A.1. The fields beyond the SM, along with their quantum numbers, are collected in table 1. The interaction between the dark sector and the SM will be given essentially from the Yukawa term (apart from the interaction term involving the standard Higgs boson and the inert Higgs doublet in the potential),

JHEP10(2016)015
acting as an inert doublet portal, where q L ∼ (2, 1/6), corresponds to the three families of SM doublets of quarks and y D is the Yukawa coupling. We will effectively consider that there is only one heavy vectorlike quark D accessible to the LHC and relevant to DM coannihilations. The possible constraints coming from these processes and also from the DM direct detection will be one of our goals. Moreover, we will choose this TeV scale heavy quark to couple only to one family of SM quarks. This choice will suppress new flavor violating effects such as on D 0 −D 0 . In particular, the case in which the exotic quark couples only to the first quark family follows by imposing minimal flavor violation: for three families of heavy quarks D iL,R with D iL ∼ d iR (D iL ∼ u iR ) and D iR ∼ q iL the spectrum for D i can be chosen hierarchical as the SM down (up) quarks and with same order and approximately diagonal Yukawa couplings (as studied in, e.g., refs. [36,37], with the difference that in our case the light-heavy quark mixing is absent due to Z D 2 ). We obtain only one heavy quark interacting predominantly to the first family after integrating out the much heavier fields. 3 The other cases are considered for phenomenological comparison.
The spectrum at the electroweak scale which we consider is an inert doublet model [12,28] augmented by an axion and a vector-like quark D interacting with the particles of the SM through eq. (2.6). The dark sector, odd by Z D 2 , consists of the fields of the inert doublet H D and the vector quark D. We choose the lightest component of H D to be lighter than D and then be part of the DM content along with the axion. It has to be noted that several models at the PQ scale can lead to this spectrum at low energies. A simple complete model that leads to this spectrum is shown in appendix A.1; it coincides with model I of ref. [15] but with a different spectrum at low energies.
The electroweak symmetry breaking is still performed by H = v/ where λ 345 ≡ λ 3 + λ 4 + λ 5 andλ 345 ≡ λ 3 + λ 4 − λ 5 . We can see that the scalar-pseudoscalar mass splitting is indeed controlled by λ 5 : In summary, the model has eight free parameters namely, where the first four elements in this set are the masses of the particles which are odd under Z D 2 , with λ 5 < 0 guaranteeing H 0 to be the lightest scalar of the dark sector besides the axion. The case in which A 0 is the lightest CP odd scalar is directly obtained replacing λ 5 → −λ 5 . As we describe in what follows, these parameters will be subjected to a multitude of constraints from the electroweak nature of the model which will reduce the viable parameter space considerably. These include theoretical constraints as well as various phenomenological ones.
Vacuum stability and perturbativity. Considerations such as vacuum stability and perturbativity restrict the range of parameters in (2.10). For the potential to be bounded from below, we need [26,38] To ensure the inert minimum ( H = v/ √ 2(0, 1) T , H D = (0, 0) T ) to be the global minimum we require [39] (scalar masses) 2 ≥ 0 , In special, the positivity of the usual Higgs mass squared requires µ 2 1 < 0. When oneloop effects are considered [40], this condition may not be strict [41]. We also require perturbativity of the scalar quartic couplings, assuming [40] |quartic self-couplings| < 4π , |X † Xhh coupling| < 4π. (2.13) Applied to the (H 0 ) 4 coupling, the first requirement 4 in (2.13) translates into λ 2 ≤ 4 3 π ≈ 4.19 [40]. A related constraint would be the unitarity in the scalar-scalar scattering matrix [42]. We do not impose the latter explicitly and argue that perturbativity already cuts off most of the non-unitary cases.
Electroweak bound. The first basic constraint comes from the electroweak nature of H D and requires that the SM gauge bosons cannot decay into the dark scalars, i.e., (2.14) LEP Limit. Susy searches at LEP [43] further exclude M H 0 < 80 GeV and LHC -dilepton + missing energy data. Using dilepton plus missing energy data from the LHC, bounds have been placed in the IDM for M H 0 < M W (the W boson mass), based on production channels such as qq In [45] the authors were able to rule out H 0 masses below 35 GeV at 95% C.L. with Run I data.
Thus far, the Higgs resonance region, where the relic density, direct, and indirect detection bounds are satisfied is left untouched. Anyway, this mass region lies outside our scenario in (2.16). (See [46] for an old study of dilepton data in the IDM).
We have reviewed the key aspects of the model as well as existing constraints for the IDM. Hereunder we discuss collider constraints based on monojet and dijet plus missing energy data from LHC at 7-8 TeV.

Collider constraints
By virtue of the Z D 2 symmetry, the vector-like quarks can only decay into a quark and a new heavy scalar, including the DM H 0 . Pair production of these new heavy quarks gives rise to SUSY-like signatures at colliders as jets plus missing energy, while associated production of a heavy quark and H 0 leads to monojets. Therefore, constraints from collider searches for supersymmetry and DM have to be taken into account prior to a dedicated study of our DM candidate. Let us discuss how we checked these collider bounds.

Bounds from SUSY and DM searches in jets plus missing energy and monojets
As aforementioned, due to the Z D 2 symmetry, the vector-like quark D can always be pair produced (DD, DD, DD) via quark or gluon fusion, or in association with a new scalar (H 0 , A 0 , H ± ) as shown in figure 1. In particular, in figure 1 we display representative contributions for pair production, diagrams (a)-(f), and single production in association with H 0 , diagram (g).
Singlet vector-quarks D can interact with the down, strange and bottom quarks via Yukawa couplings to the scalars of the model. These Yukawa couplings might be constrained by flavor physics and searches for new physics in colliders. For example, low energy physics impose constraints on the Yukawa couplings for the case where D couples with more than one family of SM quarks. We thus adopt safe benchmarks to render the model free from constraints on quark flavor violation allowing D to interact just with one family of SM quarks at a time through the Yukawa coupling y D . Figure 1. Feynman diagrams for production of DD pairs in proton-proton collisions are shown in diagrams (a)-(f). Additional diagrams obtained from crossing or charge conjugation of the initial and final state are not shown. In diagram (g) we display one contribution to D + H 0 associated production. Diagram (h) represents a subleading contribution to monojet signatures when a QCD jet from a strongly interacting line is radiated.
For the pair production of D, both QCD and Yukawa interactions with the scalars H ± , A 0 , H 0 contribute to the cross section. The t-channel diagrams with neutral scalars allow for DD and DD production alongside DD; see diagrams (d) and (e) in figure 1. A similar situation occurs in squark pair production where t-channel gluinos contribute to same-sign squarks production. Also, as in the case of squarks, the t-channel contributions impact significantly the production cross section of jets and missing energy.
It is shown in figure 2 the pair production cross section σ(pp → D 1 D 2 ) for the 8 TeV LHC for couplings with the first (down) and third (bottom) quark families, where D 1 (2) represents both a heavy quark and a heavy antiquark. The solid red (black) line represents the total cross section with all contributions from QCD and Yukawa couplings setting y D = 1 (0.5), M H 0 = 400 GeV and M A 0 = M H ± = 405 GeV. The pure QCD contribution is shown as a dashed blue line. Interestingly, the interference between the QCD and the t-channel Yukawa contributions is destructive, contrary to the SUSY case. The interference is visible only for the case of couplings with the first family, as shown in figure 2 where we can see at the left (right) panel the production cross section of vector-quark pairs with d(b)-D-H 0 coupling only. This is, of course, due to the parton content of the proton; the non-QCD t-channel diagrams connect only the initial state quarks participating in the Yukawa coupling to the vector quark D, thus, scenarios with exclusive couplings to the second and third families are suppressed and the Yukawa amplitudes contribute too little. For moderate Yukawa couplings y D 0.5, the destructive interference decreases the total cross section and only at larger Yukawa coupling regimes, where y D ≥ 1, the production rate can become larger than the pure QCD contribution.
The single and pair production of the quark D lead to monojet and two jets plus missing energy signatures at the LHC, respectively. Monojets also receive contributions from diagram (h) of figure 1 when a QCD jet is radiated from a strongly interacting particle. Monojets are striking signatures expected in the case that DM is produced in proton-proton collisions while two or more jets plus missing energy is the classical signature for production and decay of squarks and gluinos. Upper limits for production cross section times branching ratios for processes with hard jets and missing energy have been placed by the ATLAS and CMS Collaborations in the 7 and 8 TeV run of the LHC, and incorporated to the database of packages aimed to check for collider limits as SmodelS [47] and CheckMate [48].
As the quark D has a decay channel into a jet and H 0 , both the constraints from squark searches and DM searches apply in our case. In order to check these bounds we simulated the collision processes up to one extra jet to approximate higher order QCD corrections, for the 8 TeV LHC, with MadGraph5 [49], Pythia6 [50] and Delphes3 [51] after implementing the model in

JHEP10(2016)015
FeynRules [52]. Jets are clustered with the shower-k T algorithm and jet matching is performed in the MLM scheme [53] at the scale M D 4 . We checked that differential jet rate distributions are smooth across the soft-hard jet threshold.
The processes of eq. (3.1) contribute to signatures with at least two hard jets and missing energy which mimic the production and decay of squarks and gluinos. Monojet signatures receive their main contributions from the process of eq. (3.2), with a subleading contribution from eq. (3.3) where the harder jet of the event is an initial state radiation QCD jet. Experimental searches for dark matter in monojet signatures are based on exclusive criteria to select events, discarding those events with two or more harder jets [54]. For this reason, processes like eq. (3.1), with at least two hard jets, contribute little to monojets.
Collider searches constrain the parameters related to the production cross section of the process discussed above. We have chosen to constrain the Yukawa coupling y D and the vector-like quark mass M D , after fixing all the other parameters of the model. We performed scans over a wide portion of the parameters space comprising M D , M H 0 , M A 0 , M H ± and y D . For each of these points we generated 10 4 events for further analysis. The parameters scans were made as follows: We found that all scanned points passed the monojet constraints from CheckMate, but not from searches for hadronic decays of squarks and gluinos. SmodelS decomposes the full model into simplified model spectrum topologies taking into account efficiency selection criteria in order to make the correct comparison with its internal database. After that, it seeks for an experimental bound on the cross-section times branching ratio, σ(pp → D 1 D 2 ) × BR(D 1(2) → q + H 0 ) in our case, from a list of experimental publications and conference notes. Upper limits from those experimental studies on the cross sections, σ 95% , at 95% confidence level (CL), are then compared to the simulated σ(pp → D 1 D 2 ) × BR(D 1(2) → q + H 0 ). A model is considered excluded with CL above 95%, for one or more analysis, whenever we have σ(pp → D 1 D 2 ) × BR(D 1(2) → q + H 0 ) > σ 95% , or, in terms of , if the output is r 1. The first observation that we can draw from figure 3 is that the most restrictive scenario occurs when D interacts with the third family. In this case, D has a typical SUSY signature matching with searches for direct production of bottom squark pairs which translates to harder constraints in our case. Second, the bounds for the second and third families are very weakly dependent on the Yukawa coupling, an effect that we have anticipated previously. On the hand, for Yukawa couplings between 0.2 and 0.8, approximately, first family scenarios have smaller r ratios by virtue of the destructive interference between QCD and Yukawa contributions. Although the effect is not so pronounced, for y D 0.8 we see a clear trend towards the exclusion region as, in this regime, the production cross section increases.
We also see that, in general, as the mass of D increases the production cross section drops fast as shown in figure 2, but the cut efficiency somewhat compensates for the signal decrease up to 700 GeV approximately as the jets becomes harder. For masses larger than ∼ 700 GeV, the production cross section is too low and the model evades the collider constraints unless the Yukawa coupling is larger than 1.
In the next section we present the results of our analysis of the DM candidate of the model taking into account all the collider constraints we obtained.

WIMP relic density
The abundance of the WIMP is obtained in the usual way, by solving the Boltzmann equation, with the help of micrOMEGAS [67,68]. In this realization, the WIMP is in thermal equilibrium with SM particles, i.e., the annihilation and production interactions occur at similar rates in the early Universe. Although, as the Universe expands and the temperature drops below the DM mass, they can no longer be produced, and are simply able to pairannihilate. Eventually, the expansion rate equals the rate for pair annihilation and then freeze-out is established. Thus, the larger the annihilation cross section the fewer DM particles were left-over after the freeze-out. From then on, the abundance of left-over DM particle is kept basically constant. This is the standard picture, where no coannihilations are present. For the IDM, this is not the case and coannihilations play a dominant role in the WIMP abundance.
In the IDM the H 0 pair annihilation into SM particles is of the order of 6×10 −26 cm 3 /s, for 500 GeV < M 0 H < 3 TeV [69], which would naively produce an abundance below the correct value. Nevertheless, the other inert scalars H ± , A 0 interact at similar rates with H 0 and SM particles, which makes them freeze-out at a similar time. Since they are not stable, after the freeze-out they decay into H 0 increasing its abundance to match the correct value. This mechanism was explained in detail recently in [69,70]. Thus, coannihilations are an important ingredient in the IDM in order to have a viable WIMP. The setup remains identical in the WIMP+axion framework that we will advocate, as long as the coannihilations involving the exotic quark D are suppressed (to be considered in section 4.4) and axions have an insignificant relic density.
The IDM DM phenomenology can be wisely split into three mass regimes [12,28,35,71] Low mass: M H 0 < M W . In this mass range the model resembles the singlet scalar Higgs portal DM where the main annihilation modes are into light fermions, mainly bb quarks, with annihilations controlled by the quartic coupling that mix the SM Higgs and H 0 . There, A 0 and H ± have to decouple in order to avoid direct and indirect DM WIMP searches. In summary, one needs to live at the Higgs resonance to be able to reproduce the right relic density while avoiding existing constraints [94].
Heavy dark matter: M H 0 > 500 GeV. This mass region is viable and consistent with direct, indirect and collider searches. It can reproduce the right relic density thanks to coannihilations effects involving the inert scalars as we mentioned earlier, with a mass splitting 100 KeV < M H 0 −M H + ,A 0 < 15 GeV [35]. For recent studies in this mass region we

JHEP10(2016)015
refer to [69,70]. Interestingly, almost the entire parameter space of the model is expected to be probed by the Cherenkov Telescope Array [69,70].
Intermediate mass: M W < M H 0 < 500 GeV. This mass region has been entirely excluded in the light of recent direct detection limits and relic density constraints [14]. Here, the annihilation rate into gauge bosons is very efficient leading to a dwindled relic density. It is in this precise mass region which the two component DM scenario we are advocating is most relevant. Since the WIMP share its duties with the axion, the constraints are relaxed and the total relic density of Ω total = 0.1 can be achieved, motivating our work. We will explicitly show further how this is realized.

Axion relic density
As for the axion, the key question turns out to be, when is the Peccei-Quinn symmetry broken: before or after the inflation period? If it is broken before the end of inflation, the only process relevant for axion production is coherent oscillation due to the vacuum realignment and the axion relic density is given by [11,31], where θ is the initial axion misalignment angle. Note that, if θ is of order of unity, the axion can reproduce the total relic density, Ω a h 2 ∼ Ωh 2 , only for f a ∼ 10 12 GeV. We will set θ = 1 throughout. In summary, the total relic density is given by where Ω H 0 h 2 is the relic density due to the WIMP, and Ω a h 2 the one corresponding to the axion, which depends on the cosmological model. That said, it is a good timing to discuss the two component DM abundance in more quantitative terms.

Mixed WIMP-axion dark matter in the IDM
In order to take into account both axion and WIMP contributions to the total observed relic density 5 we have scanned the free parameter space in the range shown in table 2, always enforcing M A 0 − M H 0 100 keV to avoid the ruled out ineslatic DM regime [96,97]. The result of this scan is displayed in figure 4. There we show the relative WIMP and axion contributions to the total abundance as a function of f a . In figure 4 we have assumed that the exotic D quark couples only to one family of SM quarks at a time through y D , and concluded that the results are basically identical with a mild difference, within 3%, for the third family, as one can see in figure 2.
There important remarks are in order: (i) We can clearly see that for f a 5×10 10 GeV, we enter the WIMP dominated regime.
(ii) For 5 × 10 10 GeV f a 7 × 10 11 GeV, we have a plausible two component DM setup being able to meet Ω total h 2 = 0.12.
(iii) For f a > 7 × 10 11 GeV, we go into the axion dominated scenario.

JHEP10(2016)015
Parameter Scan range This plot visibly proves that one can successfully have a two component DM in the model. However, an important information in this two component DM scenario is the WIMP mass. That said, we display in figure 5 the fractions R X , with X = H 0 , a, of the total relic density as a function of the Peccei-Quinn scale f a explicitly showing the DM mass encoded in the curves. The fraction of relic abundance is defined as We have imposed M D > 300 GeV and the misalignment angle θ = 1. In addition, we have also considered the constraints (2.16) discussed in the end of section 2 and the restrictions showed in figure 3. The curve starting at R X > 80% represents the inert scalar H 0 abundance, while the curve starting at R X < 20% reflects the axion's. We enforced the total relic density to be Ωh 2 ∼ 0.1199 ± 0.0027 [98] throughout. We see clearly in figure 5 that the WIMP dominated regime favors heavier masses (M H 0 > 400 GeV), whereas the axion dominated one prefers M H 0 < 280 GeV. The reason why the WIMP dominated region prefers heavier masses is just a consequence of the IDM nature of the WIMP, since it is well known that for M W < M H 0 < 500 GeV the WIMP cannot produce Ωh 2 ∼ 0.1199 ± 0.0027. As aforementioned, this is no longer problematic in the light of our two component DM where the axion abundance makes up for the deficit, depending on the value of f a . In figure 4 the WIMP can account for 100% of the relic density as f a drops well below 10 10 GeV, because there we entered in the mass region M H 0 > 500 GeV where the relic density constraint is satisfied. The heavy quarks also play a role in setting the WIMP abundance through coannihilation processes, when M D ∼ M H 0 , as we will investigate in detail further.

New coannihilations with vector-like quarks
The DM phenomenology of the IDM from Peccei-Quinn symmetry differs from the IDM in two fundamental ways: (i) the presence of coannihilations involving the heavy vector-like quarks (D); (ii) the axion now contributes to the total relic density.
The new coannihilation processes involving the initial states H 0 D, A 0 D, H ± D and DD, will appear mediated by the Yukawa coupling  Figure 4. Contributions to the total relic density Ωh 2 ∼ 0.1199 ± 0.0027 [98] as a function of the PQ scale f a . The plot is the similar for our scenario in relation to the one presented in ref. [15]. We have assumed M D > 300 GeV, θ = 1, and the restrictions in eq. (2.16). The reason H 0 can meet the correct abundance is due to coannihilations involving the heavy vector-like quark.  Figure 5. Relative contribution of the inert scalar H 0 and axion to the total relic density, defined as R X , as a function of f a . The curve starting at R X > 80% represents the inert scalar H 0 abundance, while the curve starting at R X < 20% reflects the axion. We enforced the total relic density to be Ωh 2 ∼ 0.1199 ± 0.0027 [98] throughout. We have assumed M D > 300 GeV, θ = 1, and the restrictions in eq. (2.16).

JHEP10(2016)015
coupling y D . If the mass difference is sufficiently large or the Yukawa coupling is dwindled, the H 0 phenomenology remains identical to the IDM. To quantify the impact of these new coannihilation processes on the WIMP relic density of the IDM from Peccei-Quinn symmetry, we have used the scan over the free parameters showed in table 2. We have found that the coannihilation processes with the exotic quark D are negligible when M D 1.2 M H 0 and y D 0.7, so that we recover the DM phenomenology of the IDM in such a case, even though the coannihilation process DD → gg has pure gauge contributions independently of the Yukawa y D .
Generally speaking, coannihilation processes such as these only play a role if the mass splitting between the WIMP and the other odd particles is within 10-15%, due to the Boltzmann suppression, which is the reason for negligible coannihilation processes when We display in figure 6 [98].
Note that the coannihilations with the exotic quark decrease the WIMP population and increase the allowed DM mass compatible with the data. That is because the inclination of the relic density curve of H 0 depends on how efficient vector-quark coannihilations are. Thus, once we reach the overabundant regime, we can simply turn on such coannihilation by increasing y D and making the mass difference smaller, and bring down the abundance to the correct vale. In other words, we simply change the inclination of the abundance curves as can be explicitly seen in figure 6.
In particular, for y D = 1, right panel of figure 6, we can see a significant difference between the case in which ∆M = 200 GeV (red line), where the WIMP reproduce the total relic density for M H 0 ≈ 800 GeV, and the case in which ∆M = 100 GeV (green line), where the WIMP reproduce the total relic density for a larger mass of M H 0 ≈ 900 GeV. It is only for a splitting ∆M > 200 GeV that our model recovers the IDM phenomenology, where the vector-like quark coannihilations are turned off. For y D = 0.5 this mass difference is ∆M > 100 GeV. Notice that for y D = 1, the coannihilation cross sections are larger and hence a mass splitting must be mildly larger compared to the case with y D = 0.5 in order to suppress the coannihilations, where ∆M > 100 GeV suffices.
In the collider section we observed that y D > 0.8 might be problematic due to monojet and dijet plus missing energy constraints, therefore y D = 0.5 is a feasible benchmark model, where both relic density and collider constraints are satisfied as well as the direct and indirect DM detection probes addressed in the following.

Direct detection
WIMPs might also scatter off of nuclei and deposit an energy which can be measured by underground detectors such as LUX [99], CDMS [100] and PICO [101] among others [102][103][104][105][106][107], all of them using different target nuclei and readout techniques. By discriminating nuclear recoil from electron recoils, the experiments have been able to place stringent limits in the scattering cross section vs WIMP mass, capable of depositing an energy above few keV. In the IDM model, the direct detection limits from LUX, which is currently the world leading experiment, can be easily evaded by requiring the mass splitting between A 0 and H 0 to be above 100 KeV, and the coupling λ L to be suppressed with no prejudice to our reasoning. In particular, the values |λ L | 0.01 are well below the current sensitivity of LUX [94] and also the projected sensitivity of XENON1T [108,109].
In our model augmenting the IDM, we need to consider the presence of the exotic quark D which can mediate the WIMP interaction with the nucleus by s-channel and tchannel scattering with quarks/gluons, as shown in figure 1, diagrams (g) and (h). 6 Such interactions are governed by the Yukawa coupling y D and the exotic quark mass M D . When M D ∼ M H 0 , which is of interest to us since coannihilations with the D-quark become important, there is an enhancement in the cross section as a result of the inelastic regime. Taking y D 0.5, the model is consistent with the LUX bound on the spinindependent scattering cross section. Thus, from the left panel of figure 6, we can see that the model can simultaneously yield the right abundance and accommodate the LUX limit. For M D ∼ 1.2 M H 0 (∆M ∼ 0.2M H 0 ), when coannihilations are turned off, we find that for y D 0.7, the model is below LUX and future XENON1T [108] bounds. In summary, our benchmark model with y D = 0.5 is perfectly consistent with current and projected limits from direct detection.
Thus we conclude that the right panel in figure 6, where y D = 1, is excluded in the light of direct detection experiments. This conclusion shows the high degree of DM complementarity in our model. However, this holds true as long as H 0 accounts for the total DM abundance, which is not necessarily true in our model, specially when M W < M H 0 < 500 GeV. Since the direct detection limits are linearly proportional to the WIMP JHEP10(2016)015 local density, the bounds are alleviated and the model can be made compatible with direct detection in the regime where the axion makes up a large fraction of the abundance, i.e., for f a 7 × 10 11 GeV. We handpicked these two values for y D to show precisely when direct detection constraints become relevant and how our two component DM scenario plays an important role in satisfying both relic density and direct detection searches for WIMPs.

Indirect detection
WIMPs may self-annihilate producing a sizable amount of gamma-rays and cosmic-rays over the astrophysical background (see [111][112][113] for recent reviews). Searches for WIMP annihilations have been performed in several target regions such as the Galactic Center, Dwarf Galaxies, Cluster of Galaxies etc [69,[114][115][116][117][118][119][120][121][122][123][124][125]. In our model the mass of interest is hardly touched by current Fermi-LAT and H.E.S.S. limits [126], since the need for the axion to complement the WIMP under-abundance relaxes the indirect detection limits which depend on the local DM density squared.
Even assuming that H 0 makes up the entire DM of the Universe, for 500 GeV < M H 0 < 3 TeV, Fermi-LAT limits are rather weak, with H.E.S.S. ruling just a tiny fraction of the parameter space [69], unless boost factors are advocated [70]. It is worth mentioned that the Cherenkov Telescope Array might improve existing limits in more than one order of magnitude, and depending on the level of systematics uncertainties achieved [127,128], the entire model below 3 TeV might be excluded [69]. We emphasize though, that in our two component DM scenario such conclusions are strongly relaxed. In other words, our results are consistent with exclusion limits from indirect DM detection searches.

Conclusions
Since WIMPs and axions are arguably the most compelling DM candidates in the literature, we investigate the possibility of two component DM in a well motivated model, namely the Inert Doublet Model. We present a model that contains, beyond the SM fields, a scalar inert doublet, a scalar singlet hosting an axion, and a new vector-like quark D. These fields allow an implementation of the Peccei-Quinn U(1) P Q symmetry that solves the strong CP problem and gives rise to an invisible axion. The inert doublet originates a candidate for dark matter, stabilized by a natural Z D 2 symmetry remnant from the breakdown of U(1) P Q symmetry following ref. [15]. The new quark may provide a new portal to connect the SM to the dark sector, which is comprised of particles that are odd under Z D 2 transformations plus the axion.
Along with the WIMP, the new quark gives rise to signals involving jets plus missing energy, and also monojets at the LHC. In order to investigate possible restrictions on the parameter space of the model, we have studied all these potential signals at the LHC considering that the D quark couples to the WIMP and with just one of the SM families. We found that the most restrictive scenario occurs when D quark couples to the third family bottom quark. For example, such a scenario is excluded at 95% C.L for masses of the scalars H 0 , In our model, DM is composed by two components, the lightest inert scalar (H 0 ) and the axion. Within this scenario we performed an investigation on how the fractions of the DM relic abundance corresponding to the WIMP and to the axion change depending on the scale f a of the breakdown of the U(1) P Q symmetry, the mass of the WIMP, the masses of the other particles odd by the Z D 2 symmetry. For example, for values f a ≤ 10 10 GeV the WIMP would constitute essentially all the DM, with the axion being an irrelevant fraction of it. As f a increases the axion relic density raises, reaching a value equal to the WIMP relic density for f a 4 × 10 11 GeV.
In contrast with the inert Higgs doublet model, we found that in our model it is possible to have the WIMP from the inert doublet with mass in the interval 100 GeV M H 0 500 GeV, and comprising only a fraction of the total DM relic abundance. This region is phenomenologically important for direct detection experiments and LHC searches of exotic quarks and DM. In particular, we have shown that the IDM phenomenology remains unchanged when the coannihilations effects with the exotic quark are negligible and this happens for M D 1.2M H 0 . We conclude that one can have a plausible two component DM (WIMP plus axion) satisfying the relic density as well as collider, direct and indirect DM detection constraints.

Prospects
The assumption that the DM is composed by two or more type of particles impacts on the experiments searching for WIMPS and axions. For example, if the axion relic density constitute an irrelevant fraction of the DM the axion could not be direct detected in haloscopes experiments [129], but it could still be accessible in the projected experiment IAXO [130], which arises as a promising laboratory to test the model we proposed. On the WIMP side, direct future experiments with large exposure such as XENON1T [109] and LZ [131] are quite desired. Future collider constraints stemming from a possible 100 TeV collider or linear collider might also constrain the model even further [132][133][134][135][136].

A Simple UV completions
A.1 U(1) PQ breaking in the Higgs potential It is natural to expect that the U(1) PQ breaking at the high scale (larger than 10 9 GeV) induces at lower energies the operator in (2.5). We present here a simple model where that happens. To complete the model, we add another SM singlet scalar ϕ with PQ charge unity but inert (no vev). The relevant terms in the Lagrangian above the PQ scale will be We omit the coefficients for simplicity and, for definiteness, we take the exotic quark to be of charge −1/3, denoting it by D. The case of charge 2/3 is analogous. The PQ charges are given in table 3. After S acquires a vev S , the breaking U(1) PQ → Z D 2 is induced and we get effectively We assume |µ 2 ϕ | and the mass accompanying |ϕ| 2 to be much smaller than the PQ scale but much larger than the electroweak scale. The ϕ 2 term splits the complex scalar into two real scalars ϕ 1 , ϕ 2 of different masses M 1 , M 2 . Thus the terms with µ H , µ H of (A.2), which can be recast in the form (H † H D )(µ 1 ϕ 1 + iµ 2 ϕ 2 ), (A. 3) leads to the desired operator (2.5) with coefficient This model at the PQ scale is identical to the model I presented in ref. [15] which realizes U(1) PQ → Z D 2 in a KSVZ type axion model and, additionally, also generates neutrino masses radiatively. At the TeV scale, however, our focus is on a different physical spectrum where the DM candidates are the axion and the lightest neutral member of the inert doublet while the interaction of the heavy quark with the SM occurs also through the inert doublet. We should also emphasize that a different realization may lead to the same physical spectrum at the TeV scale -the SM augmented by an inert doublet, an axion and a exotic quark -but to a different particle content at the PQ scale.

A.2 Lighter exotic quark mass
For the model (A.2), it is expected that the exotic quark mass M D be at the order of the PQ breaking scale or at most few orders of magnitudes lower. To get M D at the TeV scale one has to tune the Yukawa coupling to at least 6 orders of magnitude. Here we show a variant where the exotic D quark have mass decoupled from the PQ scale and thus can be lighter.
The variant includes another exotic quark, which we keep denoting as D, while the original exotic quark is renamed as D . Thus the new exotic quark D has the same quantum numbers as D except that it is vector-like with respect to PQ symmetry: PQ(D L ) = PQ(D R ) = 1. Now D is still the quark that couples to the SM quarks but the QCD anomaly is generated by D .
The relevant Lagrangian is modified to Open Access. This article is distributed under the terms of the Creative Commons Attribution License (CC-BY 4.0), which permits any use, distribution and reproduction in any medium, provided the original author(s) and source are credited.