Direct detection of pseudo-Nambu-Goldstone dark matter in a two Higgs doublet plus singlet extension of the SM

We calculate the leading radiative corrections to the dark-matter-nucleon scattering in the pseudo-Nambu-Goldstone dark matter model augmented with a second Higgs doublet (S2HDM). In this model, the cross sections for the scattering of the dark-matter on nuclei vanishes at tree-level in the limit of zero momentum-transfer due to a U(1) symmetry. However, this symmetry is softly broken in order to give a mass to the dark-matter particle. As a consequence, non-vanishing scattering cross sections arise at the loop level. We find that the current cross-section limits from dark-matter direct-detection experiments can hardly constrain the parameter space of the S2HDM. However, the loop-corrected predictions for the scattering cross sections can be well within the reach of future direct-detection experiments. As a consequence, future phenomenological analyses of the S2HDM should take into account cross-section predictions beyond tree-level and the experimental constraints from dark-matter direct-detection experiments.


Introduction
Dark matter (DM) is still today one of the greatest mysteries in physics. Although the first hints of its existence were reported almost 100 years ago, and significant pieces of evidence have been gathered from different sources, we are nevertheless ignorant of its nature (see Ref. [1] for a review). Perhaps the simplest way to be in tune with all experimental results is to consider DM as a particle yet to be discovered. There are many ongoing experiments that can provide further directions in the search for the correct description of the DM field. However, in order to unmistakably observe a DM candidate, one needs direct detection (DD) experiments that probe the mass and couplings of the DM particle with the Standard Model (SM) particles via its interactions with known objects such as nuclei. As a DM particle interacts with nuclei, light and electric charge are emitted, providing information about energy and location of the collision. If we confine ourselves to the mass region of Weekly Interacting Massive Particles (WIMPs) the most restrictive and up-to-date constraints were obtained by the PandaX-4T [2], the XENON1T [3] and the LZ [4] collaborations. With the hope of a future DM detection, the significance of DD experiments in identifying the DM candidates points to the need of understanding in great detail the DM-nucleon cross sections in the different proposed models. On the other hand, as long as no DM particle is found at DD experiments, any new proposed model has to comply with the resulting upper limits on the DM-nucleon cross sections.
There is an interesting class of models where the SM particle content is extended to include an extra complex singlet field invariant under a softly broken global U (1) symmetry. Known as Pseudo-Nambu-Goldstone (pNG) DM models, they have the distinct feature of having a negligible DM direct detection cross section at leading order (LO) as first reported in Ref. [5]. In these models, the tree-level DM-nucleon cross section is proportional to the DM velocity which according to the experimentally gathered evidence is negligible when compared to the speed of light. Therefore, the first relevant contribution to the cross section comes from the one-loop electroweak corrections to the DM-nucleon cross section. These were calculated in Ref. [6][7][8] for this simplest version of the pNG DM model and shown to be several orders of magnitude above the tree-level result, and they were also shown to be several orders of magnitude larger than the finite-momentum-transfer corrections at treelevel [6]. In fact, the one-loop corrected scattering cross sections can be of the order of present experimental limits and even more so of future DM detection experiments in some regions of the parameter space.
One can consider many different extensions of the simplest version of the pNG DM model, that is, the complex singlet scalar extension of the SM with a softly broken U (1) symmetry. The hallmark of these models is a negligible tree-level DD cross section but it requires specific patterns of U (1) symmetry breaking. In fact, as shown in Ref. [9], for the simple complex singlet extension, soft U (1) breaking terms other than the quadratic ones may spoil the proportionality of the couplings between the pNG DM and the Higgs states to the squared Higgs masses, which in turn precludes the tree-level cancellation of the DM-nucleon cross section. However, it is straightforward to add an arbitrary number of doublets and still have a negligible tree-level cross section, provided the right pattern of soft symmetry breaking is implemented. Other, more sophisticated extensions with the exact same feature were also discussed in the literature [10,11].
In this work we consider an extension of the simplest version of the pNG DM model containing a second Higgs doublet, which entails numerous advantages over the version with only one Higgs doublet. For instance, the former has been shown to be able to realize a first-order EW phase transition [12,13] in contrast to the pNG DM model with only one Higgs doublet in which a first-order phase transition can be realized only if the scalar potential contains terms that spoil the tree-level cancellation of the DM-nucleon scattering cross sections [14,15]. The scalar sector of the model consists of two doublets and one complex scalar singlet and we will refer to this extension as S2HDM -the singlet extension of the of the two-Higgs doublet model. The DM candidate originates from the imaginary part of the singlet field while the real component of the singlet field aquires a vacuum expectaion value (vev) and mixes with the other two CP-even fields from the two doublets. While the pseudoscalar and charged sectors phenomenology is very similar to that of the 2HDM, there are now three CP-even fields and one DM particle. The model was already studied in Ref. [12,13,16]. In Ref. [12] the S2HDM was confronted with the relevant theoretical and experimental constraints, and, the interplay between the collider phenomenology and the dark-matter phenomenology was investigated in detail.
Since in the S2HDM the DD cross section is again negligible at tree-level, the leading terms are given by the one-loop electroweak contribution to the cross section. It was conjectured in Ref. [16] that current and future DD experiments will not be able to probe the S2HDM. However, our calculation of the next-to leading order (NLO) electroweak correction to the DM direct detection cross section will show that a significant portion of the parameter space will be probed in future DD experiments.
In our calculations we follow closely the procedure in Ref. [9]. The outline of the paper is as follows. Section 2 contains a brief description of the model and introduces our notation. In section 3 we calculate the electroweak corrections to the spin-independent direct detection cross section. In section 4, the results are presented and discussed. Finally, we present our conclusions in section 5.

The S2HDM
In order to define our conventions and notation, we briefly review in this section the S2HDM, a 2HDM extended by a complex gauge singlet field that is charged under a softly broken global U(1) symmetry. For further details on the model and its phenomenology we refer the reader to Ref. [12].

Model definitions and notation
The tree-level scalar potential of the two electroweak Higgs doublets φ 1 and φ 2 and the complex singlet field φ S is given by where the absence of explicit CP violation was assumed and all the Lagrangian parameters are considered to be real. Here the terms that uniquely involve the doublet fields are identical to the scalar potential of the 2HDM with a softly broken Z 2 symmetry that prevents the existence of flavour-changing neutral Higgs currents at tree-level. Depending on the assigned Z 2 charges of the fermions, this results in the typical four Yukawa types (see e.g. Ref. [17] for details). The terms involving the singlet field respect a global U(1) symmetry, which is softly broken by the term proportional to µ 2 χ , providing a non-zero mass for the pNG DM. After EW symmetry breaking the two Higgs doublets and the singlet field acquire real VEVs, about which the Higgs fields can be expanded in terms of the charged fields φ + 1,2 and the neutral CP-even (ρ 1,2,S ) and CP-odd fields (σ 1,2 and χ) as (2) The doublets vevs define the EW scale v = v 2 1 + v 2 2 ≈ 246 GeV and tan β = v 2 /v 1 . Assuming the previous vacuum configuration, the three CP-even fields mix and give rise to the three mass eigenstates h a,b,c . The mixing in the CP-even sector is described by an orthogonal transformation R such that, where −π/2 ≤ α 1 , α 2 , α 3 ≤ π/2 are the three mixing angles, and we use the short-hand notation s x = sin x, c x = cos x. In addition to the notation h a,b,c for the CP-even scalar states, we will also make use of the notation h 1,2,3 , where the mass ordering m h 1 ≤ m h 2 ≤ m h 3 is implied. Hence, the states h 1,2,3 can (in principle) correspond to any of the states h a,b,c whose decomposition in terms of the gauge eigenstates ρ 1,2,S is defined by the mixing angles α 1,2,3 as shown in Eq. (3). The charged scalar sector is left unchanged as compared to the 2HDM, including two physical charged Higgs bosons H ± with mass m H ± and two charged Goldstone bosons G ± that give rise to the longitudinal degrees of freedom of the W ± bosons after EW symmetry breaking. The pseudoscalar fields σ 1,2 produce a physical mass state A with respective mass m A , and a neutral Goldstone boson G 0 that constitutes the longitudinal polarisation of the Z boson. The imaginary component χ of the singlet field does not mix with the other scalar and pseudoscalar fields provided it does not have a vev. Note that the remaining dark CP-symmetry φ S → φ * S prevents χ from decaying, making it a suitable dark matter candidate. The previous definitions allow us to replace the Lagrangian parameters by a set of more physically meaningful parameters that we will use to sample the parameter space of the S2HDM, The relations between this set of parameters and the Lagrangian parameters can be found in Ref. [12].

Theoretical constraints
We required the tree-level scalar potential to be bounded-from-below (BfB) along all field directions by applying the boundedness from below conditions that were found in Ref. [18]. Additionally, we ensured the EW minimum to be the global minimum of the tree-level scalar potential to prevent the EW vacuum from decaying into other unphysical minima and from being potentially short-lived as compared to the age of the universe. The algorithm to verify for each parameter point whether there exists a global minimum of the potential with v 1 , v 2 , v S > 0 and vanishing charge-breaking and CP-breaking vevs is described in Ref. [12]. Finally, we required that the perturbative treatment of the model for a given parameter point is viable. Tree-level perturbative unitarity is achieved by imposing that the eigenvalues of the 2 × 2 scalar scattering matrix in the high-energy limit are below an absolute upper value given by 8π [12]. These constraints give rise to upper limits on the absolute values of the quartic couplings λ i and combinations thereof. These conditions are especially relevant when there are large mass splittings between one of the scalars h i and the heavy doublet states A, H ± .

Experimental constraints
In this section we will briefly describe the experimental constraints that we apply in our numerical discussion. For a more detailed discussion we refer to Ref. [12]. Regarding the Higgs sector, we required agreement with the 95% confident level cross-section limits from collider searches for additional scalar states by using HiggsBounds v. 5.10.2 [19][20][21][22][23][24]. Furthermore, we checked for the compatibility with the signal-rate measurements of the SM-like Higgs boson h 125 making use of the public code HiggsSignals v. 2.6.2 [25][26][27][28]. In case of DM masses below 125/2 GeV, the additional decay mode h 125 → χχ suppresses the ordinary decays of h 125 into SM final states. Constraints derived from the branching ratio of the invisible decay of the SM-like Higgs boson were indirectly imposed by HiggsSignals v. 2.6.2 through the global constraints on the measured signal rates of the SM-like Higgs. Further constraints related to the presence of additional scalar states are derived from measurements of electroweak precision observables (EWPO). We applied a two-dimensional χ 2 test to the oblique parameters S and T [29,30] and discarded parameter points for which the predicted values were not in agreement with the experimental fit result [31] at the 95% confidence level. 1 Constraints from flavour-physics observables (see, for instance, Ref. [31]) were taken into account by using a lower limit of tan β ≥ 1.5, and for the scan in type II we additionally used a lower limit on the mass of the charged Higgs boson of m H ± ≥ 650 GeV. Finally, we applied several constraints related to the presence of the dark matter candidate χ. One of the most important requirements is to ensure a DM relic abundance Ωh 2 prediction in agreement with the observations made by the Planck satellite, leading to a measurement of (Ωh 2 ) Planck = (0.119 ± 0.003) [34]. In our analysis, this value will serve as an upper limit for the dark matter relic abundance produced via the freeze-out mechanism predicted by our model, considering that a prediction lying below (Ωh 2 ) Planck would allow for other contributions (particle or astrophysical) to the dark matter relic abundance. We computed Ωh 2 using MicrOmegas [35]. Regarding DM direct detection constraints, the primary objective of this paper is to compute the DM-nucleon scattering cross section in the zero momentum-transfer at the one-loop level of perturbative expansion and to compare our results with the current limits set by XENON1T [3], PandaX-4T [2] and LUX-ZEPLIN (LZ) [4], and, in addition, with future limits projected for DARWIN [36]. Note that there are other planned direct detection experiments such as SuperCDMS [37], just to name an example. We have taken DARWIN as a prototype for future DD experiment. We finally note that we do not take into account constraints from the indirect detection of DM, because these constraints are only relevant in a narrow mass window of the DM below m χ 100 GeV [38], and the application of the indirect-detection constraints relies on a Monte-Carlo simulation which is computationally quite expensive (see Ref. [12] for details).

Calculation of DM-nucleon scattering cross section
As already mentioned above, in the S2HDM the cross sections for the scattering of the DM on nuclei vanish in the limit of zero-momentum exchange at tree-level. Our goal in this paper is to investigate whether radiative corrections to the cross section would make the DM state χ detectable at current or future DM direct-detection experiments. In this section we discuss the calculation of the radiative corrections, where we include the dominant contributions at the one-loop level stemming from diagrams with the scalar states in the loops. Our procedure is an extension of the calculation performed in Ref. [6]. In Sect. 4 we will then present the numerical discussion of the loop-corrected scattering cross sections in order to answer the question whether the presence of the pNG DM state χ is testable at DD experiments.

One-loop contributions to Wilson coefficients
The tree level diagram is just a t-channel χq → χq scattering where q is a quark belonging to the nucleon. The one-loop contributions to this process can be divided in three main contributions: upper vertex, lower vertex and mediator corrections. There are also box corrections that do not fit in this classification. Finally, although of higher order, the gluon initiated processes play a major role in the calculation. The one-loop contributions considered are the ones given by the topologies schematically shown in Fig. 1. These include only upper vertex and mediator corrections. Let us now discuss in detail why the remaining contributions were

discarded.
The tree-level χq → χq amplitude vanishes in the limit of zero momentum transfer (the explicit expression is given in App. A). Hence, the one-loop amplitude has to be finite in the same limit, that is, there is no need for a renormalization prescription nor for any counterterm. This was already proven in Ref. [6] for the singlet extension and again checked for our model. We explicitly verified the cancellation of the counterterm diagrams which can be carried out without specifying the individual counterterms, and thus in a generic fashion that is valid for all four Yukawa types of the S2HDM (counterterms insertions are shown in Fig. 2). As a consequence, the sum of all amplitudes is UV-finite (without the addition of counterterm diagrams), and the sum is also independent of the renormalization scale, which we verified numerically. Our analysis of the UV-finiteness of the one-loop amplitude is specific to the S2HDM, although we expect the same result to hold in a broad class of models which feature a vanishing tree-level amplitude in the limit of zero momentum transfer, because then there is no counterterm that could cancel a UV-divergent contribution at one-loop level.
In Fig. 3 we show the corrections on the external χ-legs. These corrections vanish in the limit of zero-momentum transfer, since the corresponding amplitudes are proportional to the tree-level amplitude which themselves vanish by means of the U(1) symmetry, such that the corresponding diagrams do not have to be considered. Finally we present in Fig. 4 the set of diagrams with all SM particles, the charged scalars and the pseudoscalar in self-energies and tadpole loops. We explicitly verified that only diagrams with the neutral CP-even Higgs bosons and the DM state χ in the loops give rise to non-zero contributions, whereas the diagrams with the fermions, the gauge bosons, the pseudoscalar, the charged Higgs bosons and their corresponding Goldstone bosons in the loop cancel due to the proportionality to the tree-level amplitude. We note that again this was also shown to be true for the complex scalar extension [6], but in the S2HDM there are new particles in the scalar sector and the proportionality to the tree-level amplitude is not obvious. As a consequence of this result, the one-loop corrections to the scattering cross section considered in our analysis are independent of the gauge fixing, which we also explicitly verified by calculating the amplitudes in the R ξ -gauge and varying the gauge-fixing parameter.
A set of diagrams that we did not take into account are the box contributions for the process χq → χq. This is not because the amplitudes are proportional to the tree-level amplitude but rather because their contribution is at least one order of magnitude smaller than the vertex and mediator contributions. This was checked for two different models [8,39,40] and is mainly related to the fact that the amplitude is proportional to product of two Yukawa couplings to light quarks.
With all the above considerations the set of diagrams that actually contribute to the one-loop cross-section is the one with the topologies depicted in Fig. 5, containing the upper h i χχ-vertex and the h i -propagators corrections. As discussed before, the only particles that have to be considered in the loops are the neutral CP-even Higgs bosons h 1,2,3 and the DM particle χ, since the diagrams with the other particles cancel each other out as a result of the U(1) symmetry.
Moreover, the above considerations lead us to consider only the effective scalar operator for the computation of the scattering cross sections of the DM on nucleons, where m q is mass of the quark, and C s q is the Wilson coefficient that is determined order by order in perturbation theory from the matching to the full model. Since one has to consider the scattering on both up-type quarks and down-type quarks, there are important differences between the different Yukawa types of the S2HDM. In the type I and the type LS (Lepton Specific), only the doublet field φ 2 is coupled to the quarks, independently of the quark flavour. As a result of the fact that the dependence on the mass of the different quarks is factored out of the Wilson coefficients C s q as shown in Eq. (5), in these types C s q is identical for all six quark flavours, i.e.
In contrast, in the Yukawa types II and F (Flipped) the doublet field φ 2 is coupled to up-type quarks, and the field φ 1 is coupled to down-type quarks. This gives rise to the fact that the amplitudes are different depending on whether the DM particle χ scatters on up-type quarks or down-type quarks. 2 Consequently, one finds two different Wilson coefficients which we denote in the following.
C II,F q here, but instead discuss their numerical impact in terms of the DM-nucleon scattering cross sections, as discussed in the following. However, we make the obtained expressions for the Wilson coefficients available to the public as Fortran and python routines. 3

From amplitudes to cross sections
The cross section for the scattering of χ on nucleons as a function of C s q can be expressed as [56] with the DM mass m χ and the nucleon mass m N , and with N = n, p for neutrons or protons. C s q are the Wilson coefficients of the effective DM-quark scattering operators as defined in Eq. (5). The factors f N T q are defined by the nucleon matrix elements as [57][58][59] representing the contributions of the quarks to the nucleon mass. Their numerical values have been extracted from lattice simulations and from data-driven methods to be [58,[60][61][62] The heavy quark contributions in Eq. (8) are determined by making use of the QCD trace anomaly that relates the heavy quark Q = b, c, t operators with the gluon field strength tensor [59] The quantity f N T g is then defined by the matrix element f N T g can be expressed in terms of the contributions of the light quarks, such that [56] The first sum in Eq. (8) running over the light quark flavours q = u, d, s contains the contributions from the scattering of χ directly on the light quarks. Here the contributions from the heavy quark flavours can be neglected because their contributions to the nucleon mass are tiny at the energy scales relevant for the DM scattering on nuclei. As described above, the heavy quark contributions will be included as gluon initiated processes making use of the QCD trace anomaly. This gives rise to the second sum in Eq. (8), which contains the contributions from the scattering on the gluons, where we take into account at leading order only the quark-mediated contributions. As a consequence of this approximation, this contribution can also be expressed in terms of the effective operator shown in Eq. (5). Here the important contributions arise from the heavy quarks whose couplings to the Higgs bosons are not suppressed by tiny Yukawa couplings. Thus, the second sum runs only over the heavy quark flavours b, c and t.
As previously discussed, in our computation of C s q we only include the numerically dominant corrections to the upper vertex h i χχ and the h i -propagator corrections, according to the strategy also applied in Ref. [6] for the pNG DM model with a single Higgs doublet. In this approximation, in type II and type F the amplitudes C s q are different for the up-type quarks q = u, c, t and the down-type quarks q = d, s, b, whereas in type I and type LS they are independent of the quark flavour. In the latter case, one can simplify Eq. (8) and write it as where f N is the nucleon form factor that was used in Ref. [6].

Numerical impact in light of current and future experiments
In this section we will present the numerical analysis of the DM direct-detection cross sections at the approximate one-loop level. We will start our discussion in Sect. 4.1 by analyzing whether our expressions for the one-loop contributions fulfil some theoretical requirements that can be derived from symmetry arguments in order to cross check our results. In the second step, we present the results of two parameter scan projections in the type I and the type II of the S2HDM with the goal of determining whether the DM scattering cross sections are sufficiently enhanced at the loop level such that the presence of the DM state χ could be probed at DM direct-detection experiments.

General considerations
Due to the large number of diagrams that give rise to finite contributions to the DM-nuclei scattering cross sections in the limit of zero-momentum transfer, as discussed in Sect. 3, the complete expressions for the loop corrections are rather lengthy and complicated, such that they can only be evaluated numerically. Nevertheless, the expressions have to fulfil some basic requirements that can be derived by means of symmetry arguments (see Ref. [6] for a discussion in the pNG DM model with one Higgs doublet). We will discuss here if these requirements are met by our result. This will also provide us with a first insight about the order of magnitudes of the cross sections that can be achieved in the S2HDM beyond tree-level. A more complete assessment of the phenomenological impact can be found in Sect. 4.2, where we will discuss two parameter scan projections in which we take into account the whole list of theoretical and experimental constraints mentioned in Sect. 2.2 and Sect. 2.3, respectively. The presence of non-vanishing corrections to the scattering cross sections at the loop level is related to the fact that the U(1) symmetry, under which the singlet field φ S is charged, is softly broken in order to give rise to a mass for the DM state χ. If the U(1) symmetry would be exact, the cancellation mechanism for the t-channel Higgs-boson exchange between χ and the quarks would hold at all orders in perturbation theory. A condition that the one-loop corrections have to fulfil is therefore that in the limit of m χ → 0, i.e. in the limit in which the U(1) symmetry is restored, the corrections have to vanish as well. On the other hand, if the DM mass becomes much larger than the masses of the Higgs bosons, i.e. m χ m h i , the cross sections become smaller as a result of the factor 1/m 2 χ in Eq. (8). In Fig. 6 we show the predictions for the cross sections of the scattering of χ on protons σ χp as a function of m χ in the type II S2HDM. We show σ χp for different values of the singlet vev v S , where the value of the latter is indicated by the color coding of the lines. The values of the remaining free parameters are given next to the plot on the right-hand side. The parameter values were chosen such that the theoretical constraints discussed in Sect. 2.2, in particular the perturbative-unitarity constraints, are respected. However, we did not apply the experimental constraints on the Higgs sector and the DM sector. Also shown with the dashed lines are the exclusion limits at the 95% confidence level from the XENON1T experiment [3] (blue), the PandaX-4T experiment [2] (red) and the LZ experiment [4], respectively, and the dotted line indicates the future projected exclusion limits from the Darwin experiment [36]. The gray shaded area indicates the neutrino floor [63]. As expected based on the discussion above, the cross sections vanish in the limit m χ → 0 independently of the value of v S . σ χp reaches the maximum value for DM masses that are close to the masses of the CP-even Higgs bosons h i . For DM masses that are much larger the cross sections drop again until they fall below the neutrino floor at m χ 10 TeV for the smallest values of v S considered, whereas for the largest values of v S the predictions are always within the neutrino floor. One can generically observe that overall larger values of σ χp can be achieved for smaller values of v S . This is due to the fact that for fixed values of the masses m h i smaller values of v S give rise to larger values of the quartic couplings λ 6,7,8 . These couplings act as the portal couplings between the visible and the dark sector, such that larger values of λ 6,7,8 give rise to larger values of the scattering cross sections. However, larger values of the quartic couplings also yield larger values of the annihilation cross sections and, therefore, smaller values of the predicted relic abundance. As a consequence, the parameter points with the largest values of σ χp can be expected to predict a relic abundance which is smaller than the measured DM relic abundance. The impact of the predicted DM density on the prospects of probing the S2HDM at DD experiments will be discussed in more detail in Sect. 4.2.
By comparing the theoretical predictions with the upper limit from XENON1T PandaX-4T and LZ, one can see that only for the smallest value of v S = 100 GeV considered here the current DD experiments have the potential of probing the S2HDM parameter space. It should be noted that even smaller values of v S , for which σ χp would become even larger, are excluded in this scenario as a consequence of the tree-level perturbative unitarity constraints. This emphasizes the importance of taking into account such theoretical constraints in order to give an accurate estimate of the maximum values of σ χN that can be achieved in the S2HDM. While the current upper limits from XENON1T PandaX-4T, and LZ barely constrain the parameter points shown in Fig. 6, large parts of the interval of DM masses that are shown can be probed in the future by Darwin. For instance, assuming a value of v S = 200 GeV, the expected limits from Darwin would exclude the DM mass range 30 GeV m χ 1.8 TeV. In general it is interesting to note that the scattering cross sections peak for DM masses of the order of the masses of the Higgs boson. The presence of the BSM Higgs bosons can be tested at the LHC if they are not too heavy to be produced. In this case the S2HDM can be probed in a complementary way by DD experiments and colliders.
Another theoretical requirement which has to be fulfilled by the one-loop corrections that we take into account is that the cancellation mechanism only holds in the limit of vanishing momentum transfer. As a result, the cancellation mechanism breaks down if the mass of one of the Higgs bosons is not much larger than the momentum that is transferred in the scattering process. In order to demonstrate that our result also complies with this condition, we show in Fig. 7 the predictions for σ χp as a function of the mass of one of the Higgs bosons h b , with the remaining parameters fixed to the values shown on the righthand side of the plot, and where we show here the predictions of the type I S2HDM. One can see that, as expected, σ χp increases drastically in the limit m h b → 0, independently of the value of v S as indicated by the color coding of the lines. As before, we applied only the theoretical constraints in order to produce the results shown in Fig. 7, whereas  the experimental constraints were not applied. This is important to note because values of m h b 125/2 GeV would be excluded due to constraints from the signal-rate measurements of h 125 in combination with the condition of not overclosing the universe [12]. As a result, although the direct-detection cross sections can be very large if the singlet-like Higgs boson (here h b ) is much lighter than 125 GeV, DM direct-detection experiments cannot provide additional exclusion limits in this region of the parameter space.
In the opposite limit with m h b m χ , one can observe in Fig. 7 that σ χp instead increases with increasing value of m h b . This behaviour has its origin in the fact that for fixed values of v S the absolute values of the quartic couplings λ 6,7,8 grow with increasing value of the singlet-like Higgs-boson mass m h b . As already mentioned, larger absolute values of the quartic couplings give then rise to larger scattering cross sections. The absolute values of the quartic couplings are ultimately bounded from above by the constraints from perturbative unitarity. The predictions in Fig. 7 are shown for each value of v S up the maximum value of m h b for which the parameter points were still in agreement with these bounds. Consequently, the maximum values of σ χp that are achieved here in a parameter region that is potentially not yet excluded by other experimental constraints are of the order of σ χp ∼ 10 −48 cm 2 , which is well within the range that can be tested at future DD experiments like Darwin.
Another interesting feature that can be observed in Fig. 7 is the appearance of blind-spots at certain values of m h b where σ χp drops to zero. Such blind-spot regions were also observed in the simpler case of the pNG DM model with only one Higgs doublet [6]. The presence of the blind-spots is a result of a cancellation between the amplitudes of different loop diagrams, giving rise to the fact that the sum of all amplitudes, and thus the Wilson coefficients C s q , vanish. For the blind-spot on the right-hand side it is easy to see that it appears at the point at which all CP-even Higgs bosons are mass degenerate, with m h a,b,c = 125 GeV. Even though it is questionable whether such a situation is phenomenologically viable in light of constraints from the LHC measurements, it is still an interesting observation that approximately massdegenerate scalar states could yield a highly suppressed DM-nucleon scattering cross section. A second blind-spot can be observed at roughly m h b ∼ 30 GeV, where the precise location depends on the value of v S . In addition, the location of this additional blind-spot also depends in a non-trivial way on the choice of the masses m h a,b,c and the mixing angles α 1,2,3 . For both blind-spots, it might be interesting to compute corrections beyond the one-loop level in order to analyze whether they would remain, in which case their presence would be related to an accidental symmetry, or whether the higher-order corrections eliminate the blind-spots, in which case their presence relies on a purely accidental choice of parameters.
In addition to the blind-spots that appear due to vanishing scattering amplitudes between the DM state χ and the quarks, as discussed above, in the type II and the type F S2HDM further blind-spots can appear as a result of a cancellation between the different terms in the sum over the quark contributions as shown in Eq. (8). As discussed in Sect. 3.1, in type I and type LS (at the one-loop level) there is only a single Wilson coefficient C I,LS q that enters in this sum. However, in type II and type F there are two independent coefficients C II,F u and C II,F d (see Eq. (7) and the related discussion) for the scattering on up-type and down-type quarks, respectively. If these two coefficients have the opposite sign, the sum in Eq. (8) can be suppressed even though the individual terms are unsuppressed.
In order to demonstrate this feature, we show in the left plot of Fig. 8 the cross sections for the scattering of the DM state χ on protons and neutrons in type I (orange line) and type II (blue lines) for a representative benchmark scenario. As before, we applied here only the theoretical constraints in order to ensure that the scalar potential is well behaved. One can see that at values of m h b ∼ 100 GeV the scattering cross sections in type II decrease by two orders of magnitude, whereas the cross sections in the type I remains almost constant. Moreover, it should be noted that in this interval of m h b the cross sections in type II are substantially different for the scattering on protons (solid blue line) and neutrons (dashed blue line). On the other hand, in type I both cross sections are practically equal, and consequently only one line for both the scattering on protons and on neutrons is shown. As a phenomenological consequence, one can notice that since different nuclei are composed out of a different number of neutrons and protons, a hypothetical measurement of the scattering cross sections on different kinds of nuclei could be utilized to distinguish between a DM candidate χ as predicted by the types I/LS or the types II,F, respectively.
The suppression of the cross sections in type II can be understood by the fact that one of the Wilson coefficients C II,F u or C II,F d changes the sign at the corresponding mass interval of h b . In the right plot of Fig. 8 we show the Wilson coefficients as a function of m h b for the same benchmark scenario as was used in the left plot of Fig. 8. As expected, one can see that one of the coefficients (C II,F d , dashed line) becomes negative in the mass range 50 GeV m h b 200 GeV, where the mass range coincides with the one in the left plot in which the cross sections in type II are strongly suppressed. Since in type I there is only one Wilson coefficient C I,LS q , which is identical to the coefficient C II,F u in type II (solid line), the change of the sign of C II,F d has no impact on the cross sections in type I. Finally, we note that the precise location of the blind-spot visible for type II and also the amount of the suppression of the cross sections depend on the nucleon form factors f N Tq , which are only known approximately as they are determined from lattice simulations and from experimental data. As a consequence, in the parameter regions in which the scattering cross sections are suppressed due to the accidental cancellation of contributions from different quark types with opposite sign, the relative uncertainty of the cross-section predictions associated to the uncertainty of the form factors should be regarded as larger compared to other parameter space regions in which no such cancellation takes place.
As a summary of the discussion in this section, one can conclude that the one-loop corrections included in our computation fulfil the theoretical requirements that can be derived from symmetry arguments, which serves as a non-trivial cross-check of our results. Moreover, we have demonstrated that the cross sections as predicted at the one-loop level can be well within the reach of future DM direct-detection experiments. It should be noted that we did not apply here the experimental constraints on the model parameters as introduced in Sect. 2.3. In order to verify whether the future sensitivity of DM direct-detection experiments is capable of probing parameter space regions that are not yet excluded by other experimental constraints on the Higgs sector and the DM sector of the S2HDM, we will discuss in the following section two parameter scans in the type I and the type II S2HDM in which the experimental constraints will be taken into account.

Parameter scans in type I and type II
In order to estimate the relevance of the loop-corrected predictions for the cross sections of the scattering of the DM state χ on nuclei, we present here the predictions in two parameter scan projections in the S2HDM type I and type II in which we take into account all the theoretical and experimental constraints discussed in Sect. 2.2 and Sect. 2.3, respectively. We note here that the Yukawa sectors of type I and type LS as well as the Yukawa sectors of type II and type F only differ in the couplings of the Higgs bosons to leptons. Consequently, the cross-section predictions for the DM-nucleon scattering in the type I are identical to the predictions in the type LS, and the predictions in the type II are identical to the ones in the type F. Accordingly, apart from the different collider constraints that have to be applied, the results using type I and II presented in the following also provide a good understanding of the importance of future DM DD experiments in the type LS and the type F.
In our scans we used values for the free parameters as shown in Tab. 1. We fixed m ha = 125.09 GeV in order to account for a scalar state that could, depending on its couplings, behave in agreement with the experimental measurements with regards to the discovered Higgs boson. The masses of the BSM scalars were scanned up to values of 1 TeV,  corresponding to a range that is potentially in reach of the LHC. It should be noted here that for the scan in type II we used a lower limit of m H ± > 650 GeV in order to bypass constraints from flavour-physics observables, whereas in type I we used a lower limit of m H ± > 150 GeV since the flavour constraints are much weaker (see also the discussion in Sect. 2.3). In combination with the theoretical constraints on the quartic scalar couplings and constraints from the EWPO, also the lower limits on the mass scale M and the masses m h b and m A of one of the CP-even scalars h b and the pseudoscalar A, respectively, are pushed to larger values in type II in order to account for the fact that the differences between these parameters and m H ± cannot be too large. The mixing angles were scanned over all physically distinguishable parameter space, and the lower limit on tan β was chosen according to constraints from flavour physics. Finally, the singlet vev v S is varied within the scan range of the BSM scalars. We note that due to its pNG nature the DM state χ can be light even though the global U(1) symmetry breaking has its origin at energy scales much larger than the TeV scale, such that also values of v S 1 TeV would be physically reasonable. However, as we demonstrated in Sect. 4.1, sizable values of the cross sections for the scattering of the χ on nuclei are present only if v S is of the order of the masses of the CP-even Higgs bosons or smaller. Therefore, for the purpose of determining the largest scattering cross sections that can be realized in the S2HDM it is sufficient to scan only a range in which v S is of the order of m h a,b,c (or smaller).
We have generated parameter points by scanning uniformly over the given parameter ranges. For each parameter point generated in this way, we have applied the theoretical and experimental constraints discussed above, and we have discarded the parameter points for which one of the constraints was violated. For the remaining parameter points, we have calculated the predictions for the DM-nucleon scattering cross sections. We have then compared the theoretical predictions against the current and future DM direct-detection constraints from the XENON1T, PandaX-4T, LZ and the Darwin experiment, respectively. Finally, we have also taken into account the predicted value of the DM density as obtained by assuming the standard freeze-out mechanism in order to answer the question whether the parameter points that could be probed by DD experiments would also predict a sizeable fraction of the measured DM relic abundance. Moreover, in case the predicted relic abundance is substantially smaller, we address how much this reduces the prospects of probing the corresponding S2HDM parameter space by means of DD experiments.
In the top row of Fig. 9 we show the scan points in type I (left) and type II (right) with the DM mass m χ on the horizontal axis and the DM-proton scattering cross section σ χp on the vertical axis. The color coding of the points indicates the value of m h S /v S , where h S is defined as the CP-even scalar h i with the largest singlet admixture given by R 2 i3 (see Sect. 2.1). Also indicated are the cross section limits at the 95% confidence level from the XENON1T [3], the PandaX-4T [2] and the LZ [4] experiments with blue, red and green dashed lines, respectively, and the projected future limits from the Darwin experiment [36] with the black dashed line. Finally, the gray solid line indicates the neutrino floor [63]. One can see that we find points which predict values of σ χp that are within the reach of Darwin, whereas the current experimental sensitivity by XENON1T, PandaX-4T and LZ are not sufficient to probe the S2HDM parameter space in a significant way. On the other hand, the largest fractions of parameter points feature values of σ χp that are substantially below the Darwin sensitivity, and many points are within the neutrino floor in which case a possible DM detection is not very promising even in the distant future. We note here that the range of the vertical axis for σ χp was set to 10 −52 cm 2 for a better visibility of the relevant range of σ χp for which there is experimental sensitivity, although there are parameter points featuring values of σ χp that are orders of magnitude smaller. Finally, we emphasize that overall larger values of the DM-proton scattering cross section are correlated with larger values of the ratio m h S /v S , which is in agreement with the observations discussed in Sect. 4.1.
The cross-section limits from XENON1T, PandaX-4T, LZ and Darwin, as shown in the top row of Fig. 9, were derived under the assumption that the DM particle under consideration accounts for the entire measured DM relic density as measured by Planck. However, in the S2HDM one can predict the DM relic abundance composed of the state χ assuming the standard freeze-out scenario, and in the parameter scans we only demanded that the predicted DM relic abundance is not larger than the measured value, thus leaving room for additional sources that contribute to the DM relic abundance. If the predicted abundance of the DM state χ is smaller than the measured value, the prospects for the DD of DM decrease, since the number of scattering events in the detector is smaller compared to the number of scattering events expected based on the measured DM density. In order to account for the impact of the predicted relic abundance, it is illustrative to compare the upper limits on the scattering cross section from DD experiments against the predicted scattering cross section σ χp times a scaling factor where (h 2 Ω) FO is the theoretical prediction for today's DM relic abundance based on the freeze-out mechanism (obtained with the help of MicrOmegas), and (h 2 Ω) Planck = (0.119 ± 0.003) is the value as measured by the Planck satellite [34].
In the bottom row of Fig. 9 we show the rescaled cross sections ξ FO Planck σ χp in dependence of the DM mass m χ for type I on the left and for type II on the right, respectively. Here the color coding indicates the value of the predicted DM relic abundance (h 2 Ω) FO . Since we demanded (h 2 Ω) FO ≤ (h 2 Ω) Planck (see Sect. 2.3), the parameter points all feature ξ FO Planck ≤ 1. Thus, compared to the plots in the upper row of Fig. 9, the points move towards the neutrino floor and away from the experimental upper limits on the DM-nucleon scattering cross section. Nevertheless, we find a mass interval 60 GeV m χ 300 GeV in which Darwin has the potential to probe the S2HDM parameter space in both type I and type II. One should note that many of the parameter points in this interval of m χ predict a sizable fraction (or all) of the measured DM relic abundance. Hence, the DD constraints will have the potential to probe regions of the parameter space that are especially interesting in view of the predictions for (h 2 Ω) FO . 4 For larger DM masses, additional DM annihilation channels, for instance into pairs of on-shell vector bosons, top quarks or Higgs bosons h i , become kinematically open.
As a consequence, in the range 300 GeV m χ 500 GeV we find a strong suppression of (h 2 Ω) FO , and therefore ξ FO Planck 1. This gives rise to the fact that in this range of m χ almost no points are found above the projected upper limit of Darwin. For values of m χ 500 GeV, one can see that parameter points featuring sizable values of (h 2 Ω) FO can be found above the neutrino floor, however also here the projected sensitivity of Darwin is small and limited to parameter points for which the DM state χ does not account for the whole DM relic abundance. Finally, we note that no large differences between both Yukawa types can be found. Accordingly, the prospects for probing the S2HDM parameter space at future DM direct-detection experiments can be expected to be fairly similar.

Conclusions
In this work we have calculated the one-loop electroweak corrections to the DM-nucleon scattering cross section in the extension of the SM with one extra scalar doublet and one extra complex scalar singlet, dubbed S2HDM. The model provides a DM candidate through a U (1) symmetry softly broken by dimension two terms. The tree-level amplitude of the DMnucleon process is proportional to the DM velocity and thus negligible for direct-detection experiments. This is a feature of a class of models with a pseudo Nambu-Goldstone boson as the DM candidate.
The calculation was carried out by two independent calculations and exact agreement was found. Furthermore, from the theoretical point of view, the Nambu-Goldstone nature of the DM particle would have to be reflected on a zero cross section in the limit where the exact U (1) symmetry is recovered. In fact, by making use of the non-linear formulation, it becomes clear that for low-energy processes, the Nambu-Goldstone nature of the DM particle is manifest due to the proportionality between its mass and the U (1) soft breaking parameter. We explicitly checked that our results vanish in the limit of zero DM mass and hence comply with the above mentioned symmetry arguments. We furthermore verified that there was no need to introduce counterterms as the process is zero at tree-level in the limit of zero DM velocity. The fact that our results agree with the features mentioned above was an excellent and non-trivial cross-check of our computations.
A scan of the model parameters has been performed taking into account all theoretical constraints that ensure that the electroweak vacuum is stable and that the perturbative treatment of the parameter points is valid. We have also taken into account the restrictions on the parameter space imposed by the most relevant experimental constraints related to the Higgs sector and the DM sector of the S2HDM. No parameter points have been found that could be probed by present direct detection experiments such as XENON1T, PandaX-4T or LZ, while at the same time predicting a sizable fraction of the measured DM relic abundance. However, we have found such parameter points within the reach of future experiments such as Darwin. For these parameter points, we have demonstrated that they are characterized by overall larger values of the ratio m h S /v S , where m h S is the mass of the Higgs boson with the largest singlet admixture, and v S is the U(1)-breaking vacuum expectation value of the radial component of the singlet field. Hence, although the direct detection relies on a purely loop-induced process in the S2HDM, the cross sections can be far from being too small to be probed in the future, in particular if there is no larger hierarchy between the electroweak scale and the energy scale at which the breaking of the global U(1) symmetry has its origin.
In our parameter scans in type I and type II we have not observed any major differences regarding the parameter space that could be probed in future direct-detection experiments. Moreover, no major difference between the maximum values of the DM-nucleon scattering cross sections that can be achieved have been found. However, contrary to the complex singlet extensions of the SM as well as vector DM models, the S2HDM has a very interesting feature -the matrix elements responsible for the scattering cross section may depend on the nucleon type. This is a consequence of the existence of four different types of Yukawa interactions, and that in the type I and LS both up-and down-type quarks are coupled to the same Higgs doublet, whereas in type II and F the down-type quarks are coupled to φ 1 and the up-type quarks to φ 2 , respectively. This feature would potentially make it possible to distinguish between the different Yukawa types of the S2HDM if in the future the scattering of DM will have been measured on different kind of nuclei.
Finally, we note that there are many other extensions with pNG bosons as dark matter candidates. The results obtained for this model show that most probably other extension with more fields will have more freedom and may lead to even larger cross sections once radiative corrections are considered.
respectively. Hence, in type I and type LS a = 2 for q = {u, d, c, s, t, b}, whereas in type II and type F a = 2 for q = {u, c, t} and a = 1 for q = {d, s, b}. Rewriting the amplitude M in terms of the squared masses, or vice-versa replacing the squared masses in terms of the Lagrangian parameters and the vevs, and by making use of the orthogonality of R, it is easy to show that M vanishes in the limit of zero-momentum exchange, i.e. t → 0.