Mixed WIMP-axion dark matter

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 $\mathbb{Z}_2$ 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.


I. INTRODUCTION
The evidence of cold dark matter (CDM) is overwhelming from the cosmological data, even though its detection and identification continue to be one of the most interesting and challenging problems today [1]. Many particle dark matter (DM) models have been proposed over the last few decades, one of the oldest of them being the weakly interacting massive particle (WIMP) model [2][3][4][5] (for reviews, see [6][7][8]). In the WIMP scenario, the dark matter relic abundance is obtained through the annihilation of dark matter particles in the early Universe with weak scale cross sections and electroweak scale masses [2,[9][10][11]. The fact that one gets new physics at the electroweak scale for a WIMP mass ∼ 100 GeV makes this scenario a very appealing solution to the dark matter problem [12].
The absence of CP violation in the strong sector of the Standard Model (SM) is another long-standing puzzle in the particle physics community [13]. The null results of the neutron electric dipole moment measurement experiments so far restrict the value of the coefficient θ QCD of the parity-violating E · B operator to be less than 10 −10 [14]. In the present form of the SM, this is a fine-tuning problem since there is no symmetry that protects such a small number from large higher-order corrections [15]. Therefore, a natural explanation of the smallness of strong CP violation is sought, and an elegant solution to this puzzle is given by the introduction of a global U (1) Peccei-Quinn (PQ) symmetry [16][17][18][19][20]. This symmetry is spontaneously broken at a scale much larger than the electroweak scale by a scalar field, with the axion as the corresponding massless Nambu-Goldstone boson of this U (1) PQ symmetry. The coefficient θ QCD is dynamic in this model and its small value is naturally attained in this way and is inversely proportional to the PQ scale. In this context, a large number of axion models have been proposed in the literature. The early PQ model, first proposed in [16] and further developed in [17,18], augments the SM with an additional complex scalar, charged under the elec- * suman.chatterjee@tifr.res.in † anirbandas@theory.tifr.res.in ‡ tousiksamui@hri.res.in § manibrata@berkeley.edu troweak (EW) symmetry. The Lagrangian is additionally invariant under a global U (1) symmetry, which is spontaneously broken at the EW scale. However, this model predicts large axion couplings and hence is ruled out by experiments [21]. To circumvent the experimental bounds, invisible axion models were proposed independently: the Kim-Shifman-Vainshtein-Zakharov (KSVZ) model [19,20] and the Dine-Fischler-Srednicki-Zhitnitsky (DFSZ) [22,23] model. The KSVZ model introduces heavy, colored, EW singlet quarks, in addition to the PQ scalar. In this model, the axions have no direct tree-level couplings with the SM fields and can only have induced couplings to the SM leptons. On the other hand, the DFSZ model introduces an additional Higgs field, along with the PQ scalar. This allows the axions to have natural couplings to leptons at tree level. Other than this, axionlike Goldstone particles arise in a multitude of other theoretical scenarios as well, e.g., majorons [24,25], familons [26][27][28], axions from string theory [29][30][31], axions from accidental symmetry breaking [32][33][34], etc. The axion field gains a small mass inversely proportional to the U (1) PQ -breaking scale, after the QCD condensation at a temperature of about T 200 MeV. In the early Universe, the axion can be produced nonrelativistically 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 [35,36]. 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 [37]. This makes the axion a very good CDM candidate [38][39][40], although the same feeble couplings make direct detection of these axions challenging [39].
Since current DM detection experiments have shown null results, one needs to look for other alternatives to the simple WIMP scenario. One such alternative is a multicomponent DM, where one component acts as a WIMP, whereas the other components might have very different interactions. These models are less constrained due to the fact that the fraction of WIMPs in these multicomponent DM has not been determined experimentally and hence is a free parameter. This has been shown to have interesting phenomenological consequences [41][42][43][44][45][46][47][48].
In this work, we study a two-component DM model consisting of a WIMP and the axion as the DM candidates. arXiv:1810.09471v2 [hep-ph] 27 Feb 2020 This type of model gives a unifying scenario where the PQ field, which is motivated to solve the strong-CP problem, and the WIMP, which is a natural solution to DM puzzle, can be accommodated in a single go [46,49]. Furthermore, these models can also be extended to include neutrino masses [49]. Hence, these models, and their variations thereof, can account for three of the most important puzzles in the SM. Possible UV completions are considered in [46]. As a simple realization of this, one can consider the KSVZ model of axion with an additional scalar field charged under the U (1) PQ [49]. 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 [50]. 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 WIMP-axion DM model, both from DM search experiments as well as collider searches. Using the latest limit on DM-nucleon scattering cross section from XENON1T×1 yr experiment data [51], 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×1 yr 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 DM-nucleon scattering cross section. As a result, while minimal scalar DM models are mostly ruled out by direct detection bounds [52], 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 [53][54][55].
The paper is organized as follows. Section II discusses the model and the different parameters involved. Section III talks about the different experimental bounds, and how they constrain the parameters of the model. In Sec. IV, we summarize the main results, and finally in Sec. V, we conclude.

II. THE MODEL
We consider the KSVZ model of the axion, where electroweak singlet quarks Q L and 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 vectorlike and hence do not introduce any chiral anomaly [56,57]. 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 [49]. The axion a is the Nambu-Goldstone mode of the scalar field ζ, which can couple to the vectorlike quarks as well as χ. As in the original KSVZ model, the axion can act as a CDM candidate [39]. The charges and quantum numbers of the new particles are listed in Table I.
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 righthanded down-type quarks in the SM. After electroweak symmetry breaking via the Higgs vacuum expectation value (VEV) v H , one has |H| = (h 0 + v H )/ √ 2, where h 0 is the Higgs boson. 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, also known 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 [58].
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 [59,60] and W and Z bosons, respectively [61]. 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. (5). In fact, if we take , from Eq. (5) it is evident that λ ζH has to be zero; i.e., the SM Higgs does not mix with ζ, as considered in [49]. 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 between the masses of the physical states 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. (1) has a residual Z 2 symmetry which stabilizes χ 1 . It should also be noted that in Eq. (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 [49,[62][63][64][65]. 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. Note that a small χ does not necessarily mean a small mass split between the two DM states. For example, we shall consider values around F a ∼ 10 10 GeV, which means the mass split is ∆ χ F a χ ∼ 1 TeV. This is much larger than the mass of the lighter state.
At this point, it is important to note that in this setup the complex scalar χ can, in principle, develop a VEV before the PQ field ζ. This can be prevented if parameters are tuned such that the VEV of χ, v χ = − µ 2 χ /2λ χ , remains smaller than that of ζ. It is always possible to tune χ and µ χ such that the mass of χ 1 remains fixed, and v χ remains below F a . This is because the only place where µ χ and F a appear is in the expression for the masses of the real scalars in Eq. (7). Since the mass difference between χ 1 and χ 2 is proportional to χ , this has the effect of changing the mass of the heavier scalar χ 2 . However, in our analysis, χ 2 is heavy enough to be decoupled from the particle content. With this setup, one can ensure that at a very high scale, around F a , breaking of the PQ symmetry occurs before the symmetry breaking of χ. After the PQ-symmetry breaking, χ 1 and χ 2 may or may not remain tachyonic depending on the choice of parameters. However, we have to ensure that none of them can develop a VEV earlier than the electroweak scale in order to make our mechanism work. Once the Higgs develops a VEV, both χ 1 and χ 2 become nontachyonic, irrespective of whether they were tachyonic or not prior to the electroweak symmetry breaking. The details of this mechanism is worked out in the Appendix.
The mass of the axion is generated through nonpertur-bative QCD effects and is inversely proportional to F a [66]; 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 [58], 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 vectorlike quarks obtain their mass m Q = f Q F a / √ 2 after ζ 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 Although sin θ is small, the first term cannot be ignored due to the large scale F a . Using the approximation for sin θ in Eq. (4), we obtain Note that in the presence of nonzero λ ζH , the hχ 1 χ 1 coupling vanishes at This is where we differ from [49], where the authors had set λ ζH = 0, which led to the vanishing of the hχ 1 χ 1 coupling at λ χH = 0 . This shift will play a crucial role in the following analysis. Using Eq. (5), λ ζ can be written in terms of m h , λ ζH , and λ H . This gives a family of solutions, satisfying Eq. (11). In Fig. 1, we show four contours of λ H in the λ ζH − λ χH plane for a given value of λ ζχ = 0.1. Any point on these hyperbolas satisfies Eq. (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 in the figure. One can, in principle, probe other values of λ H in this parameter space, but 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. (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.
(2), if λ ζ ∼ 10 −14 , then both the scalars h and σ can have a mass ∼ O(100) GeV. 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. (9) and (10)] to values 1, which makes the whole problem highly nonperturbative. 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.

III. EXPERIMENTAL PROBES OF DARK MATTER
Naturally, this model will have vast implications for dark matter search experiments. In addition, the LHC search for heavy vectorlike particles, as well as missing energy searches, will also test this model. Using FeynRules [67,68] to implement the model, we constrain it with the latest results from these experiments. Broadly, three avenues are explored: 1. DM relic density constraint, direct and indirect detection experiments set limits on the parameters connecting the dark sector with the visible sector.
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 1 Note that the results given in Eqs. (5)-(10) were obtained in the limit Fa v H . This approximation breaks down when the λs are set to such small values. Hence one has to start from the mass matrix in Eq. (2) and proceed without any approximation to arrive at this conclusion. observed scalar measured in the collider experiments from that of SM Higgs. This will also constrain the parameters of the model.

Since the masses of the DM and the vectorlike quarks
are lighter or near TeV range, they can potentially be produced at the LHC. Nonobservation of such particles will limit the model parameter space.
The rest of this section discusses these types of experimental constraints in detail.

A. 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 temperaturedependent potential for the axion is generated from nonperturbative QCD effects [66]. But the axion field does not start oscillating in the potential and remains frozen at its initial value until its mass becomes larger than the Hubble expansion rate H(t) =Ṙ/R, where R(t) 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 nonrelativistic speed. They contribute toward the CDM abundance today and their density is approximately given by [39,69], 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 analysis in this paper [70]. 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 as 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 χ .
The large mass split between the two states prohibits the possibility of coannihilation of χ 1 and χ 2 during DM abundance formation. As noted before, the mass split ∆ χ ∼ F a χ which, in the region of the parameter space of our interest, is much larger than m χ . For example, χ = 1 MeV and F a = 10 10 GeV imply ∆ χ = 1 TeV. During freeze-out of χ 1 , its typical kinetic energy is of order T χ ∼ m χ /20. Therefore, by this time the number density of χ 2 particles is Boltzmann suppressed relative to χ 1 , n 2 /n 1 ∼ exp (−∆ χ /T χ ) and hence is negligible. The DM relic abundance forms only through annihilations of χ 1 into SM particles.
We show the dependence of the χ 1 relic density as a function of its mass m χ in the left panel of Fig. 2. We used micrOMEGAs5.0 [71] to numerically compute Ω χ h 2 . The behavior for very small and large m χ can be understood as follows. For very small values of m χ (∼ 10 GeV), χ 1 can annihilate only into the lighter 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 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 schannel resonance from the h propagator. As m χ increases further from 62.5 GeV, the cross section falls leading to a 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 , as the colored Q L,R 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 (Ωh 2 ) obs in the left panel of Fig. 2 [1]. The overabundance region, shown as a gray shade, is disallowed. However, the underabundance region is allowed since the axion abundance Ω a h 2 can account for the rest of the relic. Therefore, the observed relic abundance We note that Ω χ is virtually independent of F a due to the v H /F a suppression in the couplings and mixing angle. Hence, F a is fixed by Eq. (13) via the Ω a h 2 term. In the right panel of Fig. 2, we show the variation of F a with m χ for three different values of θ a , for which the axion satisfies the relic constraint in the region, where the χ relic is underabundant. The result for different θ a is just a rescaling of the value of F a according to Eq. (12). The gap in between the allowed lines account for the overabundance of DM due to the χ.
We want to emphasize here that the interaction between ζ and χ fields does not affect the relic abundance of the axion and WIMP sectors. The χ ζ * χ 2 + h.c. terms in Eq.(1) introduce interactions between χ 1,2 and a, such as The interaction involving χ 2 is not important as its population is already Boltzmann suppressed. The interaction with χ 1 is ( χ /F a ) suppressed and, therefore, not relevant for relic calculation of either species.

B. Direct detection of dark matter particles
The DM direct detection experiments look for scattering between the DM particles and nuclei in the detector material. Any interaction between the DM and the SM quarks or 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.
The DM-nucleon scattering 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. (10)]. Note that this cancellation is entirely due to kinematics. There is neither any dynamical symmetry imposed to keep m χ around m h /2 nor any fine-tuning required. Since the enhancement of the cross sections is due to kinematics, the enhancement will remain stable under radiative corrections. Note that in this model, χ 1 forms only a fraction f χ of the total dark matter abundance [43,[72][73][74][75][76][77]: Therefore the DM-nucleon cross section needs to be rescaled by f χ before comparing it with the experimental

results.
Note that in the literature, another convention of rescaling the DM-nucleon scattering cross section given by f χ = Min[1, Ω χ /Ω c ] exists, which saturates f χ to unity in the overabundant region [78,79]. However, in the context of our model, we use the prescription given in Eq. (14). This is well justified, since there exists a concrete prediction for calculating the DM relic in this model. Particularly, when there is a global overdensity, we do not assume that any other unknown mechanism can account for the relic density locally. While the direct and indirect detection constraints would depend on the choice of f χ , considering them with the relic bounds does not yield any additional allowed regions. Hence, our definition does not affect the final results.
Presently, the most stringent bound on the DM-nucleon cross section is given by the XENON1T×1 yr data [51]. 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 . The rescaled cross section f χ σ χN as a function of λ χH is shown in the left panel of Fig. 3. The cross section has a dip at λ χH = 0.14 and increases as ∼ λ 2 χH for larger values as explained before. However, f χ on the other hand has a peak at the same parameter point because of inefficient relic annihilation due to vanishing of the hχ 1 χ 1 coupling. Additionally, it has an inverse relation with λ χH for larger values: f χ ∼ σ −1 ann ∼ λ −2 χH . Therefore, together the rescaled cross section f χσ χN does not have any features as shown in the left panel in Fig. 3.
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. This is shown in the right panel of Fig. 3, where we plot f χ σ χN as a function of m χ for λ χH = 0.14. The gray region shows the XENON1T constraint. In passing, we comment that the region around m χ ∼ m h is allowed by the relic constraint only if λ χH 0.079 or λ χH 0.25. However, these regions are excluded by the XENON1T bound.
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 ultralow mass. There have been a few experimental efforts to look for axionic dark matter. For example, the ADMX experiment [80] uses a rf cavity to look for its interaction with the electromagnetic field. In the KSVZ model, the interaction strength between an axion and two photons 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 a 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 [81]. 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 the correct factor to cancel out the assumption and, hence, do not require special attention.

C. Dark matter annihilation signal
Various astrophysical observations hint that the present day Universe consists of galaxies sitting inside halolike structures formed by gravitational clustering of DM particles [82]. 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%). The rescaled annihilation rate f 2 χ (σv) is shown in Fig. 4 as a solid black line. Note that here also the annihilation cross section is enhanced due to the s-channel resonance from the SM Higgs propagator. However, after rescaling with f 2 χ , which is decreased at around m χ ≈ 62.5 GeV, the annihilation rate σv shows a dip followed by a sharp increase around that point. The dependence on λ χH comes through the g hχ1χ1 coupling. The sharp dip is due to the s-channel resonance from the SM Higgs. The most stringent upper bound on this cross section is provided by the dwarf galaxy observation of the Fermi-LAT satellite data [83]. The gray shaded region is ruled out by the Fermi-LAT constraint. The green regions are excluded by the relic constraint itself, i.e., fχ > 1.
There have been many experiments which have looked for DM annihilation signals from various astrophysical objects [83][84][85][86]. 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 [83]. It is derived from 6 yr 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 Fig. 4, we show this upper bound on the annihilation cross section due to Fermi-LAT as the gray shaded region. Note that the indirect detection bounds rule out most part of our parameter space, except a region around m χ m h /2. 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. [87] within 1σ C.L. (also see [88]).

D. 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 mix 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 conjugateX, its signal strength is defined as where σ exp stands for the experimentally observed cross section of the process pp → h and BR exp is the experimentally observed branching ratio of the process h → XX. Similarly, σ SM and BR SM in Eq. (16) stand for the corresponding values predicted in the SM. We compare observed µ XX with the theoretically calculated µ XX from the model in different decay channels. Due to the mixing, the physical scalar h will have a cos θ component in all the couplings with the SM. An additional decay mode of h to χ 1 χ 1 is possible if m χ < m h /2. If the 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.  [98] In Table II, we tabulate the recent measurements of signal strength of the observed scalar h by both ATLAS and CMS Collaborations 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 Fig. 5, we show the variation of the signal strength of h in the W W * channel as a function of λ χH for two different masses of χ 1 . As expected from Eq. (17), 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. 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 (AT-LAS) 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 II 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 Fig. 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. (10)], goes to zero. This feature is evident from the plot in the right panel of Fig. 5. Away from this point, the BR increases in both sides, tending to unity for a high value of g hχ1χ1 , which indicates that Γ new is the dominant decay mode, and all other modes are suppressed. Nonobservation of this decay mode of the observed 125 GeV scalar at the LHC, therefore, places an upper limit on the invisible decays of h. These upper limits are tabulated in Table III. In the right panel of Fig. 5, the gray (green) shaded region is the area disallowed at 95% C.L. by CMS [100] (ATLAS [99]) 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 [101] and CMS [102] 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 vectorlike 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 . The presence of χ 1 will again contribute to the missing energy in the detector. The 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.
The contribution to the oblique parameters (S, T and U ) matters if the vectorlike quarks mix with our SM quarks. In that scenario, it contributes to the parameters strongly, depending on the mixing angle [103]. In our case, we do not directly have a SM top and Q L,R mixing because of the PQ symmetry. The mixing happens only through the ζ-Higgs mixing, which induces a hQ L Q R term. This mixing is suppressed by v H /F a , which is super tiny. As a result, the contributions to the S, T and U parameters are not important. 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 abovementioned final states at the 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 parton-level 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 cross sections presented in Fig. 6 have been calculated at leading order using MadGraph5 [104] with the NNPDF2.3LO parton distribution function [105]. The production cross section of QQ in the 2j+MET channel has negligible dependence on f d,s,b since the dominant partonlevel 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 generated at partonic level using MadGraph5 with the NNPDF2.3LO parton distribution function using the Universal FeynRules Object (UFO) [106] files generated by FeynRules [67,68] at a 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 [107] with 4C tune [108]. Stable particles have been clustered into anti-kT [109] jets of size 0.4 (used by both ATLAS and CMS) using the FastJet [110] software package; only the jets with P T more than 30 GeV have been considered for further analysis.
In Fig. 7, we present some important and representative differential distribution of some observables as are considered by experimental collaborations to search for signals. The top-left panel in the figure shows the distribution of p T of the leading jet while the panel in top right shows the distribution for p T of the second jet. In the bottomleft panel, we show the distribution of missing transverse energy ( / p T ). The bottom-right panel shows the distribution of H T = j ∈ jets | p Tj |, 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 Fig. 7. SM background samples have been generated with at leading order (LO) using MadGraph5 [104] with the NNPDF2.3LO parton distribution function [105] at center-of-mass energy of 13 TeV and PYTHIA 8.210 [107], with the same 4C tune [108] as used for generation of the signal sample, has been used for the simulation of fragmentation, parton shower, hadronization 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 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 [55,[111][112][113][114][115][116][117][118], 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 searches in these channels to search for dark matter signals [55,[111][112][113]. Few other models, especially supersymmetry (SUSY) in the R-parity conserving scenario, also lead to these kinds of signals. These searches have also been done by both CMS [114,115] and ATLAS [116][117][118]. 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 a search is beyond the scope of this work.

IV. RESULTS
Our main results are summarized in Fig. 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, 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 peaklike structure in Fig. 8, shows the area ruled out by DM relic abundance constraint corresponding to the 2σ bound Ωch 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 [51]. The hatched region within the red curve is ruled out by the DM annihilation data from DES-Fermi-LAT experiment [83]. 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 [89,90]. The bound coming from the signal strength of the Higgs is shown in orange [99,100]. The white, unshaded region represents the allowed parameter space in this model.
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 gammaray 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.

V. 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 a 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 vectorlike quarks at high-luminosity and high-energy operating modes of the LHC can shed further light on the viability of this model.
Higher-order loop corrections may modify the DM direct detection cross section to some extent, as discussed in Refs. [119][120][121][122]. Additionally, virtual internal bremsstrahlung in the annihilation of χ 1 may introduce special features to the gamma-ray energy spectrum, making it easier to detect by some of the experiments [123]. However, these corrections were not taken into account here and will be addressed in a future work. 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.

ACKNOWLEDGMENTS
We thank the Workshop on High Energy Physics Phenomenology 2017 for providing us with an environment for lively and fruitful discussions where this project started.
We thank Sabyasachi Chakraborty for relevant discussions in the initial stages of this project and Basudeb Dasgupta for useful inputs to the project and comments on the manuscript. We are grateful to Ranjan Laha for pointing out the effects of higher-order loop corrections in the dark matter direct detection cross section, as well as the effects of virtual internal bremsstrahlung in the observed gamma-ray spectrum. We thank the anonymous referees for the useful suggestions to improve the manuscript. We acknowledge use of the grid computing facility in Department of High Energy Physics, Tata Institute of Fundamental Research, for part of the Monte Carlo sample generation work. M. S. acknowledges support from the National Science Foundation, Grant No. PHY-1630782, and to the Heising-Simons Foundation, Grant No. 2017-228. T. S. acknowledges financial support from the Department of Atomic Energy, Government of India, for the Regional Centre for Accelerator-based Particle Physics (RECAPP), Harish-Chandra Research Institute.

Appendix: No Symmetry Breaking of χ
The potential for χ before PQ-symmetry breaking is given by The minima for this potential occurs at v χ = − µ 2 χ /2λ χ . We can always choose parameters such that v χ < F a . This will prevent χ from developing a VEV before PQ-symmetry breaking. After PQ-symmetry breaking, the potential governing the evolution of the real components of χ, viz. χ 1 and χ 2 , is given by There are four possibilities depending on the nature of the parameters: (i) µ 2 χ1 , µ 2 χ2 > 0.-This does not lead to a vev for χ 1 or χ 2 .