Pseudo-Nambu-Goldstone Dark Matter and Two Higgs Doublets

We study a dark matter model with one singlet complex scalar and two Higgs doublets. The scalar potential respects a softly broken global symmetry, which makes the imaginary part of the singlet become a pseudo-Nambu-Goldstone boson acting as a dark matter candidate. The pseudo-Nambu-Goldstone nature of the boson leads to the vanishing of its tree-level scattering amplitude off nucleons at zero momentum transfer. Therefore, although the interaction strength could be sufficient large to yield a viable relic abundance via thermal mechanism, direct detection is incapable of probing this candidate. We further investigate the constraints from Higgs measurements, relic abundance observation, and indirect detection.


I. INTRODUCTION
Astrophysical and cosmological observations suggest that the majority of matter in the present Universe consists of a nonluminous component called dark matter (DM). In the conventional paradigm, dark matter is a thermal relic remaining from the early Universe, implying that the interaction strength between DM and standard model (SM) particles may be comparable to the strength of weak interactions [1][2][3]. However, null signal results from recent direct detection experiments have put rather stringent constraints on the DMnucleon scattering cross section [4][5][6]. This has become a great challenge to the thermal DM paradigm.
In the last case, tree-level interactions of a pNGB are generally momentum-suppressed. As direct detection experiments essentially operate in the zero momentum transfer limit, the amplitude of pNGB dark matter scattering off nucleons vanishes at tree level [26]. Loop corrections could break the global symmetry, resulting in nonvanishing scattering. Nevertheless, a further investigation have shown that the DM-nucleon cross section at one-loop level is typically below O(10 −50 ) cm 2 , far away from the capability of current direct detection experiments [32]. Therefore, such a pNGB DM framework seems very appealing for thermal DM.
Previous studies in this framework assumed that the Higgs sector just involves one Higgs doublet as in the SM [25][26][27][28][29][30][31][32]. In this work, we would like to extend the study to two Higgs doublets [33]. A Higgs sector with two SU(2) L doublet has fairly good motivations. Firstly, two Higgs doublets are typically required for constructing realistic supersymmetric [34] and axion [35] models. Secondly, the flexible scalar mass spectrum and additional CP violation sources in two Higgs doublet models may be helpful for generating a desired baryon asymmetry of the Universe through the baryogenesis mechanism [36]. Finally, two Higgs doublets could provide an available portal to thermal dark matter with attractive phenomenological features [11, 18-21, 23, 24, 37-41].
In this paper, we consider that the scalar sector involves two SU(2) L Higgs doublets as well as a complex scalar S, which is an SM gauge singlet. Most terms in the scalar potential obey a global U(1) symmetry S → e iα S. The exception is a quadratic term that softly breaks this symmetry and gives mass to the imaginary part of S, denoted as χ. The real scalar χ is what we call pNGB dark matter. Its pNGB nature makes its scattering amplitude off nucleons vanish at tree level, evading direct detection constraints. Nonetheless, it is able to obtain an observed DM relic abundance via the thermal production mechanism. We will perform a random scan in the parameter space to investigate reasonable parameter points that satisfy current Higgs measurements at the Large Hadron Collider (LHC), observation of the DM relic abundance, and constraints from indirect detection experiments.
The paper is organized as follows. In Section II, we describe the details of the pNGB DM model with two Higgs doublets, including the scalar potential, mass eigenstates, four types of Yukawa couplings, the vanishing of the DM-nucleon scattering amplitude, and the alignment limit. In Section III, we perform a random scan in the parameter space and investigate phenomenological constraints from LHC Higgs measurements, relic abundance observation, and indirect detection. Section IV gives the conclusions and outlook. In Appendix A, we write down the scalar and gauge trilinear couplings. Appendix B gives some expressions for decay widths of the SM-like Higgs boson.

II. MODEL DETAILS
In this section, we study the model details. As explained above, we assume that the scalar sector involves two Higgs doublets and one SM gauge singlet, and there is a softly broken global U(1) symmetry leading to pNGB dark matter. The fermion content is assumed to be the same as in the SM. Analogous to generic two Higgs doublet models, there are four types of Yukawa couplings that do not induce flavor-changing neutral currents (FCNCs) at tree level. We find that these four types are all applicable to our purpose.

A. Scalar Potential
The two SU(2) L Higgs doublet fields are denoted as Φ 1 and Φ 2 , both carrying hypercharge +1/2. The complex scalar S is an SU(2) L singlet and carries no hypercharge. For simplicity, we make two common assumptions for the scalar potential. The first assumption is that CP is conserved in the scalar sector, leading to only real coefficients. The second one is that there is a Z 2 symmetry Φ 1 → −Φ 1 or Φ 2 → −Φ 2 forbidding quartic terms that are odd in either Φ 1 or Φ 2 , but such a symmetry can be softly broken by quadratic terms.
Under these assumptions, the general terms in the scalar potential constructed with Φ 1 and Φ 2 are given by [33] And we can write down the potential terms that involve S and respect a global U(1) symmetry S → e iα S, In addition, we introduce a quadratic term softly breaking the global U(1) symmetry, Note that even if m 2 S is complex, we can always make it real and positive by a phase redefinition of S. Then V 2 and V soft respect a dark CP symmetry S → S * [26,30]. The soft breaking term V soft can be justified by treating m 2 S as a spurion, arising from a more fundamental theory that does not induce other soft breaking terms involving odd powers of S [26,28]. Now the whole scalar potential is In particular regions of the parameter space, Φ 1 , Φ 2 , and S develop nonzero vacuum expectation values (VEVs) v 1 , v 2 , and v s . They can be expanded as while those for the CP-odd scalars are given by The above mass terms can be diagonalized by rotations where the rotation angle β satisfies Now G ± and G 0 are massless Nambu-Goldstone bosons eaten by the weak gauge bosons W ± and Z, while H ± and a are physical states with masses The CP-even scalars ρ 1 , ρ 2 , and s mix with each other. Their mass terms are − L mass,ρs = 1 2 where the elements of the 3 × 3 symmetric mass-squared matrix M 2 ρs are given by M 2 ρs can be diagonalized by a 3 × 3 real orthogonal matrix O : The mass eigenstates h i (i = 1, 2, 3) are then related to the interaction eigenstates by One of h i should behavior like the SM Higgs boson in order to be consistent with observation.
Below we adopt a convention with m h 1 ≤ m h 2 ≤ m h 3 .
From the covariant kinetic terms we derive the mass terms for the weak gauge bosons, where c W ≡ cos θ W with θ W denoting the Weinberg angle, and g is the SU(2) L gauge coupling.
Defining v ≡ v 2 1 + v 2 2 , the masses of W and Z bosons become just as in the SM. From the Fermi constant 22 GeV. Note that v 1 and v 2 satisfy v 1 = vc β and v 2 = vs β , where we have used the shorthand notations s β ≡ sin β and c β ≡ cos β.
The scalar and gauge trilinear couplings of the scalar mass eigenstates can be found in Appendix A.

C. Yukawa Couplings
Unlike the standard model, Yukawa couplings between the two Higgs doublets and SM fermions generally lead to tree-level FCNCs, which could cause phenomenological problems in flavor physics. This is because diagonalizing the fermion mass matrix cannot make sure that the Yukawa interactions are also diagonalized. Nevertheless, if all fermions with the same quantum numbers just couple to the one same Higgs doublet, the FCNCs will be absent at tree level [33,[42][43][44]. This can be achieved by assuming particular Z 2 symmetries for the Higgs doublets and fermions.
As a results, there are four independent types of Yukawa couplings without tree-level FCNCs, listed as follows. Type-II: Lepton-specific: Flipped:  Type-I Type-II Lepton-specific Flipped Thus, the interaction eigenstates u i and d i are related to the mass eigenstates u i and d i via As we would not discuss neutrino physics in this work, we assume the lepton sector is the same as in the SM.
After the scalars develop the VEVs, the Yukawa interactions provide mass terms to the fermions. For the mass eigenstates, the four types of Yukawa terms can be expressed in a same form, where P L and P R are the left-handed and right-handed projection operators, respectively.
The coefficients ξ f h i and ξ f a are listed in Table I.

D. Vanishing of the DM-Nucleon Scattering Amplitude
In this subsection, we verify that the tree-level amplitude of DM scattering off nucleons vanishes at zero momentum transfer. In our case, DM-nucleon scattering are induced by DMquark scattering. Therefore, we just need to prove that the DM-quark scattering amplitude vanishes in the zero momentum transfer limit.
From the U (1) symmetric potential (2), we obtain the trilinear couplings for the DM candidate χ as where the coupling coefficients for the mass eigenstates are given by At tree level, only the CP-even Higgs bosons h 1 , h 2 , and h 3 can mediate χ scattering off quarks. The Feynman diagram is shown in Fig. 1.
Take the Type-I Yukawa couplings as an example. Defining a Lorentz invariant t ≡ p µ p µ , where p µ is the 4-momentum of the mediator h i , we can write down the DM-quark scattering amplitude as where u(k 1 ) andū(k 2 ) are the wave functions for the incoming and outgoing quarks, respectively. In the zero momentum transfer limit, t → 0, and the above amplitude can be re-expressed as . From Eqs. (31) and (20), we have Utilizing these equations as well as the orthogonality of O, we obtain This can be understood as the amplitude expressed in the interaction basis [26]. The inverse of M 2 ρs can be expressed as its adjugate A divided by its determinant, i.e., (M 2 ρs ) −1 = A/ det(M 2 ρs ). The relevant elements of A are We then have Therefore, we have proven that the tree-level DM-quark amplitude iM vanishes in the zero momentum transfer limit for the Type-I Yukawa couplings. Similarly, we can prove this for the Type-II, lepton-specific, and flipped Yukawa couplings.
As the global U(1) symmetry is softly broken, loop corrections would give a nonvanishing DM-nucleon scattering cross section [26]. Nonetheless, we expect that the loop-induced cross section should be typically O(10 −50 ) cm 2 , as suggested by the one-loop evaluation in Ref. [32] where only one Higgs doublet is considered. Thus, current and near future direct detection experiments should not be able to probe our pNGB DM model.

E. Alignment Limit
Current LHC Higgs measurements favor a 125 GeV SM-like Higgs boson. If one of the CPeven Higgs bosons mimics the SM Higgs boson, the constraints from Higgs measurements can be easily satisfied. For the two Higgs doublets, such a situation can be achieved by requiring the additional scalars are much heavier than the weak scale so that the lightest CP-even Higgs boson reproduces SM-like Higgs signals at the LHC. This is known as the "decoupling" limit [45]. In general, a particular parameter set or relation leading to a CPeven Higgs boson mimicking the SM Higgs boson is referred as an "alignment" limit. The decoupling limit is of course an alignment limit, but it is less interesting, as the new particles might be too heavy to be accessed at the LHC.
A more interesting possibility is alignment without decoupling [46,47]. In order to find such a possibility, we may rotate the two Higgs doublets Φ 1 and Φ 2 into the Higgs basis [48,49] Φ and have Now Φ h gains a VEV v and contains a CP-even scalar h as well as the Nambu-Goldstone bosons, while Φ H have zero VEV and contains a CP-even scalar H and the physical states H + and a. Consequently, the tree-level interactions of the CP-even scalar h with weak gauge bosons and SM fermions are totally identical to those of the Higgs boson in the SM. Therefore, the alignment limit means that h does not mix with H and s.
In the Higgs basis, the potential terms (1) transform to where the new parameters are related to the previous parameters by [40] On the other hand, the potential terms (2) transform to where the new parameters are given bỹ Then the stationary point conditions for the scalar potential are As a result, the mass-squared matrix for CP-even scalars (h, H, s) is In order to prevent h-H and h-s mixings, the off-diagonal terms (M 2 hHs ) 12 and (M 2 hHs ) 13 should be absent, corresponding toλ This is the alignment condition in our model. When this condition is satisfied, the tree-level couplings of h to SM particles are exactly the same as those of the SM Higgs boson.

III. PHENOMENOLOGICAL CONSTRAINTS
In this section, we take the Type-I Yukawa couplings as an illuminating example to investigate the phenomenological constraints from Higgs measurements, relic abundance observation, and indirect detection.
Then we require the selected parameter points must give positive m 2 h 1,2,3 , m 2 H + , and m 2 a , ensuring physical scalar masses. Moreover, one of the CP-even Higgs bosons h i should have a mass within the 3σ range of the measured SM-like Higgs boson mass m h = 125.18 ± 0.16 GeV [50]. We recognize this scalar as the SM-like Higgs boson, and denote it as h SM , and further examine if its properties are consistent with current measurements.
In the κ framework [51], the couplings of the SM-like Higgs boson to SM particles can be expressed as where g SM hgg , g SM hγγ , and g SM Zγγ are the loop-induced effective couplings to gg, γγ, and Zγ, respectively. κ's are coupling modifiers, whose values are all equal to 1 in the SM. Eq. (A6) implies that κ W and κ Z are equal in our model, and we will use κ V representing both of them. Assuming the SM-like Higgs boson is h SM = h i , we have The coupling modifiers for fermions can be read off from Table I. For the Type-I Yukawa couplings, all SM fermions have the same coupling modifier, given by It is also helpful to define another modifier κ H as where Γ SM h is the Higgs total decay width in the SM, Γ h SM is the total decay width of the SM-like Higgs boson h SM , and Γ BSM h SM is the h SM decay width into final states beyond the SM (BSM). Thus, κ H indicates the deviation of the Higgs width decaying into SM final states and is also equal to 1 in the SM. In our model, Γ BSM h SM can be generally separated into two parts:  [52] to study the constraints from current Higgs measurements. Lilith is able to construct an approximate likelihood based on experimental results of Higgs signal strength measurements. For each selected parameter point in our random scan, we put the corresponding m h SM , κ V , κ f , BR inv , and BR und into Lilith. Then Lilith can evaluate κ g , κ γ , and κ Zγ involving the loop contributions from SM fermions and gauge bosons whose couplings are modified by κ f and κ V , including NLO QCD corrections. Such an evaluation have neglected the loop contributions from the BSM scalars in our model. Nonetheless, these scalars are typically heavy and/or have small couplings. Therefore, their contributions are insignificant for most of the selected parameter points. We further use Lilith to calculate the likelihood −2 ln L for each parameter point based on Tevatron data [53], ATLAS Run 1 data [54][55][56][57][58][59][60][61], CMS Run 1 data [62][63][64][65][66], ATLAS Run 2 data [67][68][69][70][71][72][73][74][75], and CMS Run 2 data [76][77][78][79][80][81][82]. We then transform −2 ln L to a p-value and require that the selected parameter points should give p-values larger than 0.05. This means that we have rejected the parameter points that are excluded by data at 95% confidence level (CL). Now we can analyze the properties of the remaining parameter points. Fig. 2 shows the Lilith p-values of the selected parameter points projected in the tan β-λ 1 and tan βλ 2 planes. We find that when tan β 0.2 (tan β 5), λ 1 (λ 2 ) tends to converge on λ SM = m 2 h /v 2 0.26, which is the quartic Higgs coupling in the SM. This is because tan β i.e., Φ 1 (Φ 2 ) acting as the SM-like Higgs doublet. Since experimental data favors an SM-like Higgs boson, the corresponding quartic coupling would close to its SM counterpart.
Additionally, we project the parameter points in the m h SM -m h 1 and m h 2 -m h 3 planes in Figs. 3(a) and 3(b), respectively. In Fig. 3(a), the points with h SM = h 1 align along a horizontal line with m h 1 125 GeV, while the remaining points indicate that the SM-like Higgs boson is not the lightest CP-even Higgs boson h 1 . On the other hand, two sets of aligned points in Fig. 3(b) correspond to h SM = h 2 and h SM = h 3 .
The projection on the m H + -m a plane is presented in Fig. 4. From Eq. (16), we know that the difference between the masses of the charged Higgs boson H + and the CP-odd Higgs boson a are due to the λ 4 and λ 5 couplings. Ifm   and λ 5 contributions, the difference would be negligible, as demonstrated in Fig. 4 for m H + , m a 500 GeV. Fig. 5(a) shows the projection on the |κ V |-|κ f | plane. We find that the parameter points with |κ V | |κ f | 1 have the largest p-values, implying that current data still favor that the 125 GeV Higgs boson has SM-like couplings. Nonetheless, |κ V | may range from ∼ 0.85 to ∼ 1, and |κ f | may range from ∼ 0.6 to ∼ 1.3. In addition, there are two categories of parameter points approximately aligning along two outstanding lines. and Eqs. (60) and (61) become κ V κ f O 2i , where κ V and κ f have a nearly total positive correlation. As |O 2i | ≤ 1, in this case both |κ V | and |κ f | cannot exceed 1. Most of the parameter points in this category corresponds to the horizontal line with |O 2i |/s β 1 in the |O 1i |/c β -|O 2i |/s β plane shown in Fig. 5(b), while the rest gives |O 2i |/s β < 1.
There are some scatter points not belonging in the two categories. Most of them correspond to tan β ∼ 1.
The dominant contributions to κ g come from the top and bottom loops, leading to a parametrization of [50] κ On the other hand, κ γ is mainly contributed by the W and top loops, resulting in [50] In both cases, the interference between the two contributions gives a term with negative coefficient. In Fig. 6(a), we project the parameter points in the κ g -κ γ plane, where the points also align along two lines. One line implies a positive correlation between κ g and κ γ , corresponding to Category 1. This is because the relation κ V κ f gives rise to such a positive correlation via Eqs. (64) and (65). On the other hand, when |κ V | 1, κ γ is negatively correlated to |κ f | as all selected parameter points satisfy κ V κ f > 0. As κ g is positively correlated to |κ f |, Category 2 results in a second line with a negative correlation between κ g and κ γ .
κ Zγ is also dominantly contributed by the W and top loops, given by [50] The correlations of κ Zγ to κ V and to κ f are similar to those of κ γ . In addition, κ H can be expressed as [50] κ 2 where all the coefficients are positive. Thus, κ H is positively correlated to both |κ V | and |κ f |. The projection in the κ Zγ -κ H plane are shown in Fig. 6(b). Analogous to Fig. 6(a), Category 1 leads to a line indicating a positive correlation between κ Zγ and κ H in Fig. 6(b). Besides, parameter points in Category 2 roughly align along a second line with a negative correlation.
In Fig. 7(a), we show the projection in the m χ -BR inv plane. When m χ > m h SM /2, we have BR inv = 0, because the invisible decay h SM → χχ is kinematically forbidden. When m χ < m h SM /2, the invisible branching ratio BR inv could be as large as ∼ 25% and still consistent with data at 95% CL. The projection in the Γ h SM -BR und plane are presented in Fig. 7(b). We find that the undetected BSM branching ratio BR und can be allowed up to    ∼ 27%, while the total width Γ h SM can range from ∼ 2 MeV to ∼ 7 MeV. There is a line implying a positive correlation between Γ h SM and BR und . This is reasonable, because opening new decay channels enlarges the total width.
In order to investigate the alignment limit, which corresponds toλ 6 =κ 1 = 0, the selected parameter points are projected in the tan β-λ 6 and tan β-κ 1 planes in Figs. 8(a) and 8(b), respectively. We find that most of the selected points satisfyκ 1 0, showing no particular dependence on tan β. On the other hand,λ 6 is typically close to 0 for tan β 20 and tan β 0.05. For 0.05 tan β 20, there is no particular favor in the alignment limit.

B. DM Relic Abundance
The thermal relic abundance of dark matter is essentially determined by the total velocityaveraged annihilation cross section at the freeze-out epoch, which we denote as σ ann v FO . In our model, the DM candidate χ has the following annihilation channels if kinematically allowed.
• Annihilation into a pair of fermions, χχ → ff . This channel is mediated by s-channel CP-even Higgs bosons and suppressed by fermion masses. Thus, tt and bb are the important final states.
• Annihilation into a pair of weak gauge bosons, χχ → W + W − , ZZ. This channel is also mediated by s-channel CP-even Higgs bosons.
• Annihilation into a weak gauge boson and a Higgs boson, χχ → W ± H ∓ , Za, mediated by s-channel CP-even Higgs bosons.
• Annihilation into a pair of CP-even Higgs bosons, χχ → h i h j . This channel can be mediated by s-channel CP-even Higgs bosons, as well as by t-and u-channel χ.
Additionally, there are contributions from quartic scalar couplings.
• Annihilation into a pair of CP-odd or charged Higgs bosons, χχ → aa, H + H − . This channel is contributed by the mediation of s-channel CP-even Higgs bosons and quartic scalar couplings.
Some numerical tools are adopted to calculate the relic abundance. We implement the model with FeynRules 2.3.34 [83], and import the generated model files to a Monte Carlo generator MadGraph5 aMC@NLO 2.6.5 [84]. Then we utilize a MadGraph plugin MadDM 3 [85] to compute the relic abundance Ω χ h 2 for each parameter point.
The relic abundance predicted by the selected parameter points is shown in Fig. 9, where the color bar denotes the freeze-out annihilation cross section σ ann v FO . We find that the observed value Ωh 2 = 0.1186 ± 0.0020 given by the Planck experiment [86] corresponds to σ ann v FO ∼ O(10 −26 ) cm 3 /s, which is typical for thermal dark matter. Increase in m χ typically reduces the annihilation cross section, and hence increase the relic abundance. Consequently, if the DM candidate is too heavy, say m χ 3 TeV, the observed relic abundance could not be achieved.
In Fig. 9, the parameter points predicting Ω χ h 2 over the observed value by 2σ are denoted with crosses. These points are considered to be excluded by data, because DM overproduction by the thermal mechanism contradicts standard cosmology. On the other hand, if the predicted thermal relic abundance is too low, there could be some nonthermal production occurring after DM freezes out.  [86]. The parameter points denoted with crosses lead to DM overproduction, while those denoted with dots are allowed by data.

C. Indirect Detection
In this subsection, we discuss constraints from γ-ray indirect detection experiments. There are couples of dwarf spheroidal galaxies discovered as satellites of the Milky Way Galaxy. They are considered as the largest substructures of the Galactic dark halo, predicted by the cold DM scenario [87,88]. As known so far, they are the most DM-dominated systems [89]. Moreover, γ-ray emissions from typical astrophysical sources, such as neutral and ionized gases, and recent star formation activity, are expected to be rare in such dwarf galaxies [90][91][92]. These properties make them perfect targets for searching for γ-ray emissions from DM annihilation.
The DM velocity dispersion in dwarf galaxies is typically ∼ O(10 −5 ) [93], which is smaller than DM velocities at the freeze-out epoch by four orders of magnitude. Therefore, if the velocity dependence is significant in DM annihilation, the total velocity-averaged cross section in dwarf galaxies σ ann v dwarf could be much different from the freeze-out value σ ann v FO .
We further use MadDM to calculate σ ann v dwarf for each parameter point assuming the average DM velocity is 2 × 10 −5 . The ratio of σ ann v dwarf to σ ann v FO is demonstrated in Fig. 10(a), where the parameter points excluded by the Planck relic abundance measurement are not shown. Most of the parameter points give σ ann v dwarf / σ ann v FO ∼ 2. Nonetheless, some points give the ratio away from O(1), indicating significant dependence on velocity. This is typically due to DM annihilation through the resonances of CP-even Higgs bosons, since the resonance effect extremely depends on the difference between the resonance location and the velocity-dependent center-of-mass energy [94,95].
The vertical dashed line in Fig. 10(a) indicates the location of m χ = m h SM /2, correspond-  ing to the resonance of the SM-like Higgs boson. We can see that the ratio around this line could range from ∼ 10 −4 to ∼ 30. On the other hand, locations of the other resonances are not fixed, but their effects are also important. Fig. 10(b) shows the projection of the parameter points in the m χ -σ ann v dwarf plane, as well as the 95% CL upper limits on σ ann v dwarf given by a analysis of Fermi-LAT and MAGIC γ-ray observations [96]. The analysis combined 6-year observations of 15 dwarf galaxies from the Fermi-LAT satellite experiment and 158-hour observations of a single dwarf galaxy Segue 1 from the MAGIC Cherenkov telescopes. The limits were obtained assuming that DM solely annihilates into bb. However, there are various DM annihilation channels in our model. Fortunately, the γ-ray spectra yielded from these channels should be similar to the spectrum from the bb channel, because they are contributed by similar processes, such as hadronization, hadron decays, and final state radiation. Therefore, we have a good reason to expect that the bb limits are approximately applicable to our case. From Fig. 10(b), we can observe that a large fraction of the parameter points with m χ 1 TeV are ruled out, while the parameter points with m χ 100 GeV and Ω χ h 2 ∼ 0.1 are not excluded. Additionally, if m χ m h SM /2, the resonance effect could both yield a data-allowed relic abundance and lead to a small σ ann v dwarf evading the indirect detection constraint.

IV. CONCLUSIONS AND OUTLOOK
In this paper, we have studied the pNGB DM framework with two SU(2) L Higgs doublets Φ 1 and Φ 2 . The DM candidate χ is the imaginary part of a complex scalar S, which is an SM gauge singlet. Most of the scalar potential terms respect a global U(1) symmetry S → e iα S, expect for a soft breaking term giving mass to χ. As a result, χ becomes a stable massive pNGB. Mass eigenstates in the scalar sector also include three CP-even Higgs boson h i , a CP-odd Higgs boson a, and charged Higgs bosons H ± .
There are four possible types of Yukawa couplings without tree-level FCNCs, just as in usual two Higgs doublet models. DM scattering off nucleons is mediated by the CP-even Higgs bosons. Because of the pNGB nature of χ, the scattering amplitude vanishes in the limit of zero momentum transfer for all the four Yukawa coupling types. Although loop corrections lead to a nonvanishing amplitude, the scattering cross section is expected to be below ∼ O(10 −50 ) cm 2 , Consequently, current and near future direct detection experiments are incapable of probing such a DM candidate.
Taking the Type-I Yukawa couplings as an example, we have performed a random scan in the 13-dimensional parameter space. The selected parameter points are required to provide an SM-like Higgs boson whose properties are consistent with current LHC Higgs measurements. We have found that for tan β 1 or tan β 1, one of the Higgs doublets acts as the SM-like Higgs doublet, and data favor the alignment limit. On the other hand, for tan β ∼ 1 there is no preference to the alignment limit.
We have also calculated the relic abundance and annihilation cross sections predicted by the selected parameter points. For m χ 3 TeV, it is possible to achieve the observed relic abundance. Because of the resonance effect, the present velocity-averaged annihilation cross section at dwarf galaxies could be rather different from that in the freeze-out epoch. Fermi-LAT and MAGIC observations of dwarf galaxies have excluded a large fraction of parameter points with m χ 1 TeV. Nonetheless, for m χ m h SM /2 or 100 GeV m χ 3 TeV, it is still possible to simultaneously satisfy the constraints from the relic abundance observation and indirect detection.
Such a pNGB DM model is strongly related to Higgs physics. The proposed future Higgs factories, such as CEPC [97], ILC [98], and FCC-ee [99], would greatly improve the Higgs measurements. We expect that these measurements could significantly restrict the parameter space in our model. Nevertheless, Higgs measurements are not able to pin down the DM candidate mass m χ , which is solely determined by the soft breaking term that does not affect the rest scalar masses. Thus, indirect detection experiments in the future are essentially important for exploring this model. where s W ≡ sin θ W . The derivative symbol ← → ∂ µ is defined as F ← → ∂ µ G = F ∂ µ G − G∂ µ F . The coupling coefficients are given by This appendix gives the decay widths of the SM-like Higgs boson into two-body BSM final states when they are kinematically allowed. Assuming the SM-like Higgs boson is h SM = h i , its invisible decay width at tree level is Moreover, its decay widths into aa and H + H − are given by Furthermore, the h i → aZ decay width can be expressed as where the λ function is defined by λ(x, y, z) ≡ x 2 + y 2 + z 2 − 2xy − 2xz − 2yz.
The decay width of h i → H + W − is given by which is equal to the decay width of h i → H − W + .
If h SM = h 2 or h 3 , it is possible to decay into h 1 h 1 and h 2 h 2 , whose widths can be commonly expressed as withg ijj = g ijj + g jij + g jji . If h SM = h 3 , there is another possible decay channel into h 1 h 2 .