Electroweak Corrections in a Pseudo-Nambu Goldstone Dark Matter Model Revisited

Having so far only indirect evidence for the existence of Dark Matter a plethora of experiments aims at direct detection of Dark Matter through the scattering of Dark Matter particles off atomic nuclei. For the correct interpretation and identification of the underlying nature of the Dark Matter constituents higher-order corrections to the cross section of Dark Matter-nucleon scattering are important, in particular in models where the tree-level cross section is negligibly small. In this work we revisit the electroweak corrections to the dark matter-nucleon scattering cross section in a model with a pseudo Nambu-Goldstone boson as the Dark Matter candidate. Two calculations that already exist in the literature, apply different approaches resulting in different final results for the cross section in some regions of the parameter space leading us to redo the calculation and analyse the two approaches to clarify the situation. We furthermore update the experimental constraints and examine the regions of the parameter space where the cross section is above the neutrino floor but which can only be probed in the far future.


Introduction
Ever since Dark Matter (DM) became an inevitable ingredient in model building, all kind of proposals integrating DM candidates into phenomenologically viable models have emerged, from the most simple extensions of the Standard Model (SM) to fairly intricate models. The accumulated data from astrophysics and cosmology strongly suggests that if DM is a particle, it is most probably cold and with a mass close to the electroweak scale. Particles with these features are known as Weakly Interacting Massive Particles (WIMPs). In this work we study a simple extension of the SM where a complex singlet is added to the SM field content. The model is built such that after spontaneous symmetry breaking one of the singlet components will mix with the SM-like Higgs boson while the other one will play the role of the DM candidate. This type of extension was first proposed in Refs. [1][2][3][4] and recently reviewed with updated experimental constraints in Ref. [5].
A particular version of this extension known as the Pseudo Nambu-Goldstone DM model (PNGDM) has a scalar potential invariant under a global U (1) symmetry which would give rise to a Nambu-Goldstone boson. The symmetry is then softly broken and a pseudo Nambu-Goldstone boson emerges as the Dark Matter candidate. As discussed in previous works [6,7], the nature of this particle makes the DM-nucleon tree-level cross section proportional to the velocity of the DM particle and therefore negligible (see also Ref. [8] for an interesting discussion on the subject of Goldstone and pseudo-Goldstone bosons). Hence, the leading order cross section is given by its one-loop contribution. The first calculation of the electroweak corrections to this process was performed in Ref. [7] and almost at the same time a second calculation appeared in Ref. [9]. In this work we will perform once again the calculation of the electroweak corrections to the DM-nucleon cross section while discussing in detail the main differences with respect to the two previous calculations and the reasons for settling this issue. Our calculation will be performed with a different renormalisation scheme. This allows us to perform a rough estimate of the remaining uncertainties on the cross section due to missing higher-order corrections.
We then perform a scan in the parameter space taking into account the most relevant theoretical and experimental constraints. We will show that there are still allowed points in the parameter space above the neutrino floor [10] but only experiments in the far future will be able to probe them.
The paper is organized as follows. In section 2 we briefly present the complex singlet extension of the SM while in section 3 we introduce the renormalisation of the model. The various aspects of the DM direct-detection cross section at tree level and at one-loop level are discussed in 4. In section 5 we discuss the results and future prospects of DM detection in this model taking into account the most recent constraints. We summarise our findings in section 6.

The Model
A simple extension of the SM by a scalar gauge singlet is enough to provide a valid DM candidate. The new complex scalar field S, with zero hypercharge and zero isospin only enters the model via the scalar potential that can be written as where the mass parameters µ 2 H , µ 2 S , m 2 χ and the quartic couplings λ S , λ H , λ HS are real due to hermicity. The doublet H and singlet S fields are expanded as follows with the electroweak vacuum expectation value (VEV) v and the singlet VEV v s . With this definition the model is invariant under the DM charge conjugation S → S * , which guarantees the stability of the imaginary part of S. Furthermore, in order to simplify the potential, an invariance under the Z 2 symmetry S → −S has also been imposed. The real part of S develops a vacuum expectation value (VEV), while the doublet develops the usual (SM) VEV that gives mass to the SM fermions and gauge bosons, with v = 2m W /g, m W the W boson mass and g the SU (2) coupling constant. Because the real part of S acquires a VEV, Φ S cannot be a viable DM candidate and it mixes with the doublet real neutral component Φ H . Using the minimum conditions we can write the mass matrix of the two neutral states as The mass eigenstates h 1 and h 2 are obtained from the gauge eigenstates via with the orthogonal matrix R(α) R(α) ≡ cos α sin α − sin α cos α .
One of these mass eigenstates is identified as the 125 GeV Higgs boson. The DM particle is given by χ. Exploiting the tadpole conditions Eq. (4) its mass can be written as We also require that the potential is bounded from below inducing the tree-level conditions The parameters of the potential can be written as functions of the masses, mixing angle and the VEVs as We choose the following parameters as independent input parameters, Some final comments regarding the scalar potential are in order. The potential is invariant under a U (1) symmetry (S → e iα S) that is softly broken by the dimension-two term proportional to m 2 χ . The Goldstone boson related to the U (1) symmetry acquires a mass proportional to m 2 χ . Due to the Z 2 symmetry there are no more terms contributing to the mass of the pseudo Nambu-Goldstone boson. Hence, the U (1) symmetry is recovered by setting m 2 χ = 0, where the true Nambu-Goldstone boson is recovered.

Renormalisation of the PNGDM
In the following, we present the renormalisation of the PNGDM in order to be able to calculate the electroweak (EW) corrections to the scattering process of the pseudo-Nambu Goldstone DM particle with a nucleon. Having defined the full set of input parameters in Eq. (11), the bare parameters p 0 are replaced with the renormalized ones, p, according to where δp corresponds to the counterterm of the respective bare parameter p 0 . For a generic bare field Ψ 0 (scalar, fermion or vector field), the renormalized field Ψ is expressed as with the field strength renormalisation constant Z Ψ ≡ 1 + δZ Ψ . Note that Z Ψ is a matrix for mixing fields as present in the PNGDM. Dropping for simplicity the index 0 in the following, we hence make the following replacements and analogously for the tadpole parameter In the following we discuss each sector separately.

Gauge Sector
The gauge sector of the PNGDM is not extended compared to the SM. To set our notation and conventions we will list the counterterms of the gauge sector in the following. We choose to perform the renormalisation in the mass basis of the PNGDM, so that the following set of on-shell (OS) counterterms are taken for the gauge sector e → e + δZ e e , (17c) As m W and g have already been defined above: where m Z is the mass of the EW neutral gauge boson Z and the electric coupling is denoted by e. The renormalized fields are obtained through the field strength renormalisation constants as Applying OS conditions yields the following mass counterterms with T indicating the transverse part of the self-energies Σ ii (ii = W, Z). The counterterm for the gauge coupling g, is obtained from the one for the electric charge and the one for the Weinberg angle θ W using The electric charge counterterm is fixed in the Thomson limit, which by making use of Ward identities allows us to write [11] 1 where we introduced the short-hand notation s W ≡ sin θ W , c W ≡ cos θ W , and The corresponding wave-function renormalisation constants guaranteeing the correct OS properties are given by

Scalar Sector
In the PNGDM we have two additional scalars, one extra CP-even Higgs boson and the DM candidate χ. The two CP-even scalars are mass-ordered as h 1 and h 2 with m h 1 < m h 2 and the SM-like Higgs boson, with a mass of 125.09 GeV, can be either of them. We again use an on-shell scheme for the fields. The field strength renormalisation constants read The mass matrix with the additional tadpole contributions is given by The rotation matrix R(α) is defined in Eq. (7) and diagonalises the gauge eigenstates in the Higgs mass basis. The tadpole terms δT in the tree-level mass matrix are bare parameters. At next-to-leading (NLO) they get shifted due to EW corrections to the vaccuum state of the potential. Defining the tree-level VEV to be the same to all orders of perturbation theory, requires the introduction of tadpoles counterterms such that the the one-loop renormalized onepoint functionT i (i = H, S) vanisheŝ Note that the rotation matrix from the gauge states to the Higgs mass states also applies to the tadpoles, yielding the relation between the tadpoles T i (i = H, S) and T h i (i = 1, 2) The one-loop mass counterterm of the Higgs sector is then given by with Equation (29) is strictly expanded to one-loop order, so that terms O (δαδT i ) are dropped. Applying OS conditions yields (i = 1, 2) There is just one DM candidate and therefore the renormalisation constants are defined by with the self-energy Σ χχ of the DM candidate.

Quark Sector
In the quark sector we assume a diagonal CKM matrix for simplicity. This means we neglect flavor mixing and the OS scheme is applied for each quark individually. The field strength renormalisation constant has to be formulated for the left-and right-handed field of the quarks (q = u, d, s, b, t) [12] q L/R → 1 and the mass counterterm is introduced through The two-point correlation function of the quarks is written as where the superscripts L, R and S correspond to the left-, right-handed and scalar parts of the self-energies, respectively. The ω ± are the left-and right-handed projectors. The full set of counterterms is then given in terms of the left-/right-handed and scalar parts of the respective self-energies as

Renormalisation of the Mixing Angle
The rotation Eq. (6) of the interaction states Φ H and Φ S to the mass eigenstates h 1 and h 2 introduces the mixing angle α that needs to be renormalized as well. The renormalisation of the mixing angles in SM extensions was thoroughly discussed in Refs. [13][14][15][16][17][18][19][20][21][22][23][24][25]. There are many possibilities to renormalize the mixing angle. One possibility is to use a physical process, like a decay. However, it is known that the usage of a process-dependent scheme may yield an unphysically large counterterm [14] which in turn leads to extremely large corrections. In this work we will use the KOSY scheme, proposed in Refs. [26,27], which connects the angle counterterm with the usual OS counterterms for the scalar. 2 The bare parameter α 0 can be expressed in terms of the renormalized one, α, as Considering the field strength renormalisation before the rotation, and expanding it to strict one-loop order, yields the field strength renormalisation matrix √ Z H connecting the bare and renormalised fields in the mass basis. Using the rotation matrix expanded at one-loop order results in Demanding that the field mixing vanishes on the mass shell is equivalent to identifying the off-diagonal elements of √ Z H with those in Eq. (24), With Eq. (32) the mixing angle counterterm reads We do not perform a comparison of various renormalisation schemes, like a process-dependent, MS, or the KOSY scheme, in this work. We note, however, that in our previous work [28], when comparing these three schemes, we found that only the KOSY scheme led to reasonable NLO predictions.

Renormalisation of the Singlet VEV v s
In the Standard Tadpole (ST) scheme that we are using in this work, there is no need to renormalize the singlet VEV v s . It was shown in Ref. [29] that when choosing an R ξ gauge in the ST scheme there is no divergence associated with v s at one-loop order if the scalar field obeys a rigid invariance. In these SM extensions the singlet field is disconnected from the gauge sector and hence invariant under global gauge transformations. This is exactly the case for typical extended scalar sectors with a singlet field, like the complex (or real) singlet extension of the SM or the Next-to-2-Higgs-Doublet Model, where a real singlet field is added to the 2-Higgs-Doublet Model. We note, however, that in the alternative tadpole scheme as defined in Ref. [30] for the SM and in Ref. [17] for the N2HDM this is no longer true and a counterterm for v s is needed.

Spin-Independent Cross Section in the PNGDM
In the following, the calculation of the spin-independent (SI) cross section for the direct detection of DM is presented. The starting point is the scattering process of a DM particle with the nucleon. The effective coupling of this process is denoted by α n , where it is additionally assumed that the momentum of the nucleon is not altered, that is, the momentum transfer between the DM particle and the nucleon is negligible. We can then use the normalisation for the spinors, u n u n = 2m n . The DM-nucleon cross section with the interaction in Eq. (45) is then given by where m n corresponds to the nucleon mass and m χ to the DM mass. Since the nucleon is a bound state the contributions to the effective DM-nucleon coupling is on the one hand given by the light valence quarks (q = u, d, s) and on the other hand by the gluon interactions. In order to calculate the cross section the parton basis is used to describe the interaction between the DM and the nucleon. The SI DM-nucleon cross section is calculated by taking the related operators in the non-relativistic limit. The parton operator basis forming the most general SI-interactions for scalar DM is given by [31] with the operators The operators are built with the DM field χ, the quark spinor q and the gluon field strength tensor G µν a and α s denotes the strong coupling constant. The operator O q S describes the interaction induced by the quark-DM interactions and O g S the one induced by the gluon-DM interactions. The twist-2 operator O q µν also contributes to the SI cross section due to additional gluon induced interactions. Assuming on-shell nucleon states |n , the expectation values of the operators in Eq. (48) can be expressed as [32][33][34] n| m qq q |n ≡ m n f n q , with the nucleon matrix elements f n q and f n g calculated on the lattice. The numerical values for the matrix elements used in the analysis are given in Sec. 5.1. Eq. (49) allows to formulate the effective DM-nucleon coupling α n in terms of the Wilson coefficients defined in Eq. (47) as [35] The numerical values of the second momenta of the quarks q n (2) are also given in Sec. 5.1.
In order to give an estimate of the DM-nucleon cross section the remaining task is to calculate the Wilson coefficients C q/g S and C q T in Eq. (50).

SI Cross Section at Tree Level
We will start by showing that the SI cross section vanishes in the limit of vanishing momentum transfer. The Feynman diagrams representing the quark contributions together with the corresponding amplitude are given by where m h i (i = 1, 2) are the neutral Higgs boson masses and t = (p χ − p q ) the Mandelstam variable. The momenta of the DM particle and the quark are denoted by p χ and p q , respectively. The amplitude in Eq. (51) allows to read off the Wilson coefficient C q S by identifying m qū u as the operator O q S . The amplitude and therefore the Wilson coefficient is proportional to the momentum transfer t and vanishes in the limit of vanishing momentum transfer. Hence, the quark contribution to the SI cross section is zero. Note that this behaviour is related to the U (1) symmetry of the model as will be discussed later. Let us show that also the gluon part of the cross section vanishes in the same limit which implies that the SI cross section vanishes at tree level in the limit of vanishing momentum transfer. (Note that the twist-2 operator does not contribute at leading order.) The QCD trace anomaly allows to relate the quark operator of the heavy quarks Q = b, c, t with the gluon field strength tensor yielding the effective gluon interaction with DM particles The corresponding Feynman diagram is depicted in Fig. 1 and can be calculated by first calculating the process with a (heavy) external quark as in Eq. (51) and using Eq. (52) to determine the effective gluon interaction. These amplitudes are then used to the extract the Wilson coefficients C g S . Note that the gluon contributions are extracted in the same way as in Eq. (51), hence proportional to the momentum transfer t. Consequently, also the gluon contributions vanish in the limit of vanishing momentum transfer.

EW Corrections to the SI Cross Section
As shown in the last section, the SI cross section of the PNGDM is suppressed at tree level due to its proportionality to the momentum exchange. Since we work in the limit t = 0, the tree-level cross section vanishes and we have to calculate the cross section in the next order of perturbation theory. In this section we calculate the EW corrections to the SI cross section which follows very closely our approach presented in Ref. [28] and updated in Ref. [36]. Note that the vector DM model presented in our previous work does not show the tree-level suppression present in the PNGDM, hence the cross section is now calculated by taking the NLO amplitude squared whereas in vector DM model the LO times NLO term was taken. The generic one-loop EW

Mediator Corrections
In this section we will discuss the propagator corrections. To calculate the one-loop corrections to the mediator we first evaluate all genuine one-loop diagrams in Fig. 3 and construct the corresponding counterterm. This can be achieved by evaluating the renormalized one-loop propagator (i, j = 1, 2) with the renormalised self-energy matrix We now have everything to determine the contribution of the diagrams in Fig. 3. Note that the field strength renormalisation constant δZ is introduced artificially, since the Higgs bosons correspond to an internal degree of freedom. As it turns out, if the field strength renormalisation constants are included in all separate topologies (lower vertex, upper vertex and mediator corrections), they cancel each other exactly in the sum. Hence, in the end no artificially introduced δZ parts remain in the calculation. The inclusion of these δZ factors on the other hand allows to check for the UV finiteness in each topology by itself, simplifying the calculation or rather the bookkeeping of the contributions.

Upper Vertex (upV) Corrections
The upper vertex corrections -referred to as upV -are depicted in Fig. 4. Diagrams a) to f) are the genuine one-loop corrections and are calculated in the limit of vanishing momentum hj Figure 3: The one-loop EW corrections to the mediator. They can be split in the genuine one-loop diagrams diagrams a)-e)) and the respective counterterm amplitude (diagram f)). The indices i, j, k, l = 1, 2 indicates the respective Higgs mediator h1, h2. The possible field insertions are given by Φ1 = {hi, χ, G (0,±) , Z, W ± } and Φ2 = {χ, G (0,±) , Z, W ± , ηZ , ηW , f }, where f stands for all SM fermions, G (0,±) for the neutral and charged Goldstone bosons, respectively, and ηZ,W for the ghost fields.
transfer (i.e. incoming momentum is equal to the outgoing momentum, p in = p out ). Note that this specific limit is stricter than taking q 2 = (p χ − p q ) 2 = 0 implying for instance vanishing Gram determinants complicating the reduction to the standard one-loop scalar integrals. The numerical evaluation of the integrals is performed with the Collier package [37][38][39] and explicitly cross-checked with an in-house implementation. The counterterm diagram Fig. 4 g) is constructed by varying the tree-level coupling of the χχh i vertex (i = 1, 2) yielding Note that, since we are using the standard tadpole scheme, the introduction of a counterterm for the singlet VEV v s is not required to obtain a UV finite result (cf. discussion in 3.5). The corresponding counterterm amplitudes read then with the quark Higgs coupling (i = 1, 2) and the quark spinors u. The artificially introduced δZ factors for the internal Higgs mediator are again included to ensure the proper cancellation in the sum of all topologies.

Lower Vertex (loV) Corrections
In  The quark q corresponds to the up-or down-type quark depending on the field insertion, respectively. Note that for simplicity a diagonal CKM matrix is assumed.
tree-level coupling of the lower vertex qqh i is given in Eq. (58), hence the counterterm for this vertex reads and the full CT amplitude The presence of charged particles in the final states indicates additional infrared (IR) divergencies in the amplitudes. The introduction of real radiation to regulate these IR divergencies does not work in this context, since the matching to the parton operators in Eq. (48) occurs at the amplitude level and the cancellation of the IR divergencies happens at the cross section level. Furthermore, the inclusion of real corrections would also introduce additional tensor structures in the amplitude which have to be accounted for in the parton operator basis. The IR divergent parts of the amplitude form a closed subset of diagrams referred to as QED subset in the following and all diagrams contain an internal photon line. The corresponding diagrams are depicted in Fig. 6 where the self-energy of the quarks enters through the mass counterterm δm q and the field strength renormalisation constants δZ L/R qq , and the vertex corrections are part of the genuine one-loop vertex corrections of the lower vertex. One possible solution is the expansion of the QED subset in terms of the external quark momentum p q yielding an IR safe result as discussed in Ref. [36]. However, the U (1) symmetry of the potential leads to the complete cancellation of the QED subset, such that no IR divergencies are present in the final renormalized amplitude of the lower vertex corrections. Hence no additional treatment is required to regulate IR divergencies. The box and triangle topologies contributing to the DM-quark interactions are presented in Fig. 7, where the incoming momenta are denoted by p 1 and p 2 , respectively. For simplicity, the triangle diagrams containing Goldstone bosons (G 0 , G ± ) are not shown, but they are included in the calculation treated in the same way as the Higgs mediator triangle diagrams in Fig. 7. The definition of the momenta reflecting the vanishing momentum transfer limit allows to express the diagrams as

Box Diagrams
The generic couplings are defined as A ij = a i a j b i b j and B ij = a i a j b ij , where a i,j and b i,j are the coefficients of the h i,jq q and h i,j χχ vertices, respectively, and b ij is the coefficient of the h i h j χχ vertex. The coefficients are given explicitly by The main contributions to the amplitudes in Eq. (62) come from the regions close to the poles of the propagators, that is where k 2 is close to the squared Higgs masses m 2 1 and m 2 2 which are of the order of several hundreds to thousands of GeV 2 . In direct detection experiments, the target nucleus is almost at rest and hence the energy of the nucleons can be approximated by the Fermi energy, which is in the order of MeV. Therefore the approximation p 2 k is valid in these integrals and the denominators that contain p 2 can be expanded as follows [40,41], and using the Dirac equation / pu(p) = m q u(p) we obtain The expanded amplitudes in Eq. (65) can then be reduced with standard techniques to the Passarino-Veltmann integral basis. Furthermore, we emphasise that the expansion leads to reduced scalar integrals not depending on kinematic variables as s allowing to use the matching procedure to the parton operator basis.

General Mapping to the Wilson Coefficients
All diagrams of the NLO corrections presented in Secs. 4.2.1 to 4.2.4 have only two independent spinor structures contributing to the SI cross section, namelyū(p 2 )u(p 2 ) (with the remainder of the amplitude independent of momenta) and terms proportional to (p 1 · p 2 )ū(p 2 ) / p 1 u(p 2 ). Hence, the amplitude can be cast into the following form with some momentum-independent constants A and B. The definition of the twist-2 operator allows to reformulatē where the asymmetric part does not contribute to the SI cross section and it can therefore be dropped. The resulting amplitude and the coefficients can be mapped to the effective Lagrangian containing the parton operators Identifying the coefficients in Eq. (68) with the Wilson coefficients in Eq. (47) yields By using Eq. (69) the calculated renormalised amplitude can be mapped to the corresponding Wilson coefficient allowing to determine the SI cross section at NLO.

Gluon Contributions
Besides the DM-quark interactions also the DM-gluon interactions contribute to the SI cross section. As shown in Sec. 4.1 the leading DM-gluon contributions can be obtained by using the relation between the heavy quark operators and the gluon field strength tensor in Eq. (52), but this returns a vanishing SI cross section for vanishing momentum transfer (t → 0). Therefore, next-to leading order effects have to be taken into account to determine the DM-gluon interactions. The leading non-vanishing gluon interactions are shown in Fig. 8 which are 2-loop diagrams. The first two diagrams correspond to the generic mediator and upper vertex EW corrections in combination with the effective vertex ggh i which can be calculated in the heavy quark limit (by using Eq. (52)). The third term corresponds to an effective two-loop calculation which will be discussed later. The first two diagrams can be calculated using the renormalized upper vertex (Sec. 4.2.2) and the mediator corrections (Sec. 4.2.1) with external quarks instead of the gluon and using Eq. (52) to effectively determine the DM-gluon interactions. By identifying the gluon parton operator O g s the respective Wilson coefficient can be deduced in accordance with the quark operators. This method of including the gluon contributions poses several problems, however.
The first problem is that, as will be shown latter, the correct mass dependence, in the limit m χ → 0, is not recovered in the limit of zero DM velocity. As discussed in Ref. [8] for an exact symmetry the Goldstone boson completely decouples from all of its interactions in the limit of vanishing momentum. Furthermore, it can be shown with the help of a toy model that scattering amplitudes involving Goldstone bosons vanish in the zero-momentum limit although this is not manifest at the Lagrangian level and only occurs through a nontrivial cancellation of terms in the S-matrix. The reason is that in the zero-momentum limit the Goldstone state is a symmetry transformation of the ground state and therefore indistinguishable from the vacuum in this limit [8]. The pseudo Goldstone case is similar -we just have to take simultaneously the limit of zero-momentum together with m χ → 0 which takes us back to the potential invariant under U (1). The second problem is that by this matching, only the diagrams with electroweak corrections to the Higgs boson propagator and the upper DM-Higgs boson vertex can be taken into account. However, electroweak corrections to the lower quark-Higgs boson vertex would obviously interfere with the quark triangle, which makes a matching to heavy quarks non-trivial, since the loops do not factorize. This second problem reveals itself in the framework of our renormalisation scheme in the following two points. By neglecting the lower vertex corrections, the cancellation of the artificially introduced δZ h i h j does not occur anymore and an uncanceled finite piece of δZ h i h j remains in the Wilson coefficient. In particular, the off-diagonal elements of the wave function renormalisation constants δZ h i h j (i = j) introduce a mass pole 1/ m 2 h 2 − m 2 h 1 yielding a parametric enhancement for nearly degenerate mass spectra. This divergent enhancement does not correspond to a physical phenomenon but rather to a wrong method for the determination of the DM-gluon interactions. Also, the KOSY scheme for the renormalisation of the mixing counterterm δα produces numerically stable (in the sense of no unphysical parametrically enhanced EW NLO corrections or divergencies) NLO predictions if either δα and δZ h i h j appear in a specific combination or if δα appears in a full process several times canceling the mass pole structure [21]. The former occurs e.g. in 1 → 2 Higgs decays yielding a δZ h i h j for the on-shell Higgs state and a corresponding δα counterterm in the vertex counterterm. The latter is present for instance in the 2 → 2 scattering process χq → χq, since δα comes both from the upper and from the lower vertices. By neglecting in the lower vertex the EW corrections in the triangle-type diagrams in Fig. 8 the conditions for a numerically stable KOSY mixing counterterm δα are not given, and hence a non-physical enhancement is expected. The third term in Fig. 8 corresponds to an effective two-loop calculation, where the two different gray blobs are explained in Fig. 9. The diagrams in Fig. 9(a) are calculated using the approach presented in Ref. [40] and already applied to the VDM in Ref. [28]. Applying the heavy quark limit (valid for mediator masses below the top quark mass) allows us to formulate an effective vertex h i h j gg where a i,j are the Higgs-quark couplings defined in Eq. (63). The vertex is produced by the effective Lagrangian and therefore the Wilson coefficients O g S can be extracted by calculating the one-loop diagrams induced by the vertices depicted in Fig. 9(b). The one-loop corrections induced by the last vertex in Fig. 9(b) have to be calculated with caution. The first two vertices do not yield a UV pole in the amplitude, hence no counterterm is required. On the other hand, the last vertex generates in general a UV pole requiring a vertex counterterm, since these corrections correspond to an effective vertex correction. However, the U (1) symmetry of the PNGDM ensures the cancellation of all UV poles, yielding a UV safe amplitude.
We emphasise that the inclusion of such effective vertex corrections has to be done with caution, since the cancellation of the UV poles is not guaranteed and is model dependent. Furthermore, these corrections are effective two-loop calculations, where other two-loop contributions are dropped because they are assumed to be small. This is not the case in general. Nevertheless, the size of the included effective two-loop corrections is small compared to the other EW NLO corrections (upV,loV,med,box) when a scan over the allowed parameter space is performed. Hence, we have included these corrections in our calculation. In the following we will refer to the inclusion of the EW NLO corrections of the upper vertex or mediator in combination with the effective Higgs-gluon vertex as the approach with the additional gluon contributions. Whereby, the proper SI cross section is calculated solely by taking the effective two-loop contributions into account (third diagram of Fig. 8. As we will discuss later, these contributions yield only a sub-percentage effect on the overall cross section, hence the inclusion of these contributions does not alter the results significantly.

Numerical Set-Up and Parameters
In the following we list the numerical values used for our study. The SM input parameters are taken as [42] The SU (2) electroweak gauge coupling g and the Weinberg angle are expressed in terms of the gauge boson masses and the electroweak VEV Note that we chose to renormalize the Higgs sector in the mass-ordered Higgs basis h 1 and h 2 with the masses m h 1 < m h 2 . One of the Higgs bosons is identified as the SM-like Higgs boson with a mass of [43] m h = 125.09 GeV , and the non-SM like Higgs boson will be referred to as φ, with mass m φ . Both mass hierarchies m h < m φ and m φ < m h are allowed in the analysis. In the following, we refer to the SI cross section as the SI cross section obtained by the scattering on a proton where the proton mass is given by The nuclear matrix elements for the proton needed in Eq. (50) are [32][33][34] and it should be noted that the uncertainties in the determination of these nuclear matrix elements are not taken into account. For the parameter region scan we implemented the PNGDM in ScannerS [44,45] which is now publicly available 3 .
The points generated using ScannerS have to be in agreement with the most relevant experimental and theoretical constraints. ScannerS allows to check that the potential is bounded from below, that there is a global minimum and that perturbative unitarity holds. The SM-like Higgs couplings to the remaining SM particles are all modified by the same factor. Hence, the bound on the signal strength [43] is used to constrain this parameter. There are new contributions to the massive gauge-boson self-energies, Π W W (q 2 ) and Π ZZ (q 2 ). The variables S, T, U [46,47] are used to guarantee agreement with the electroweak precision measurements at the 2σ level.
The collider bounds from LEP, Tevatron and the LHC are all encoded in HiggsBounds 5.6.0 [48] and HiggsSignals 2.3.1 [49]. Agreement at the 95% confidence level is asked using the exclusion limits for all available searches for non-standard Higgs bosons, including Higgs invisible decays. The corresponding branching ratios are calculated with AnyHdecay 1.1.0 [45]. This code includes the Higgs decay widths, including the state-of-the art higher-order QCD corrections, for the complex singlet model as obtained from sHDECAY [50]. The code sHDECAY is based on the implementation of the singlet models in HDECAY [51,52]. For our calculations all EW radiative corrections in HDECAY are turned off for consistency.
The DM relic abundance for each model is calculated with the MicrOMEGAs code [53], which is compared with the current experimental result (Ωh 2 ) obs DM = 0.1186 ± 0.002 from the Planck Collaboration [54]. We do not restrict the DM relic abundance to be exactly at the experimental value but rather that the value predicted by the model has to be equal to or smaller than the observed central value plus 2σ. This way, we can consider both the dominant and subdominant DM cases simultaneously. Regarding direct detection the XENON1T [55,56]  (78)

Results and Discussion
We start the discussion with the Xenon plot in Fig. 10. The effective SI DM-nucleon cross section is shown as a function of the DM mass m χ . Note that the actual SI cross section has to be rescaled with the factor with the observed relic density Ωh 2 DM and the produced relic density Ωh 2 χ for the DM WIMP χ. As discussed, we do not demand that the DM candidate accounts for the full relic density. When DM is under abundant, the effective factor in Eq. (79) corrects the cross section accordingly. The relic density is calculated in the standard freeze-out mechanism with the help of MicrOmegas implemented in ScannerS. The color code in Fig. 10 denotes the value of the non-SM like Higgs boson mass m φ and the gray shaded region corresponds to the neutrino floor [10]. The different lines correspond to the limits of the different DM detection experiments. The vertical red line indicates the half of the SM-like Higgs boson mass.
All parameter points shown in Fig. 10 are compatible with the theoretical and experimental constraints described previously. The figure shows that for the entire range of the DM mass from roughly 40 GeV up to 1 TeV, only small mass regions around m χ ≈ m h /2 and m χ ≈ m φ /2 may yield an effective SI cross section above the neutrino floor. In the case of m χ ≈ m h /2 we can see a large number of points that are basically above but close m h /2; points below m h /2 are excluded by the LHC invisible decays constraints. For the region where m χ ≈ m φ /2 only a few points for m φ of the order of 1 TeV are above the neutrino floor. There are however many points in this region that are above the region where most points are concentrated. The fact that only scattered points appear in this region is related to a combination of the experimental constraints. These regions correspond to the two resonances h and φ, respectively. The requirement of proper dark matter abundance leads to the suppression of the coupling between DM and the resonance. However, the kinematical enhancement caused by the resonance compensates for the suppressed couplings that govern DM annihilation in the early Universe. Parameter points below the neutrino floor are not of interest, since those points wiil not be able to be checked by future direct detection experiments, as the neutrino floor puts a natural limit to the sensitivity of this kind of experiments. The abrupt cut for m χ below m h /2 is induced by Higgs to invisible searches yielding a strict limit, since in this parameter region the decay h → χχ is kinematically allowed. Hence, only a few allowed points are found in this specific parameter region. We emphasise that the tree-level prediction for the SI cross section is zero due to the vanishing momentum-transfer limit, hence the parameter points cannot be constrained by direct detection experiments with tree-level calculations. However, as shown in Fig. 10 the EW NLO corrections can shift the parameter points above the neutrino floor and in the reach of the future Xenon 100T experiment. Therefore, the EW NLO corrections might play an important role in the discussion of the sensitivity of the direct detection experiments and derived exclusion limits. In Fig. 11 the SI cross section is shown as a function of the non-SM like Higgs boson mass m φ with the color code indicating the mixing angle α. Note that we do not include the factor f χχ here. The SI cross section drops for degenerate neutral Higgs boson masses (m φ ≈ m h ) because the NLO cross section is proportional to m 2 φ − m 2 h as shown in Ref. [7]. On the lefthand side, all of the about 260.000 parameter points fulfilling the theoretical and experimental constraints are plotted. On the right-hand side, only the parameter points that lead to direct detection cross sections above the neutrino floor are plotted in color and all remaining parameter points are shown in gray. It is interesting to note that there are allowed points with very large cross sections which, however, do not fulfil the relic density constraints. This way most points with the appropriate relic density have a cross section below ∼ 10 −46 cm 2 , except for a few very heavy non-SM like Higgs boson masses. All parameter points with an SI cross section above the neutrino floor have a maximal mixing between the Higgs doublet gauge state Φ H and the singlet Φ S . Only a single parameter point is above the neutrino floor with one neutral Higgs boson being a singlet-like Higgs boson, meaning that the mass eigenstate is almost given by the singlet field component. This parameter point is also the only parameter point having an inverted Higgs spectrum (m φ < m h ) while providing an SI cross section above the neutrino floor.
In Fig. 12 we show the SI cross section as a function of the DM mass m χ , where the DM mass is varied while keeping the other input parameters fixed. On the left side the resulting SI cross section is shown by starting from the benchmark point given in Tab. 1 and then varying only the DM mass while keeping all other parameters fixed, and on the right side we show the results by starting from the benchmark scenario presented in Ref. [9] with the input parameters and variable DM mass. The green line corresponds to the SI cross section calculated in the approach presented in Sec. 4.2 and the blue line shows the result for the approach with the additional inclusion of the gluon contributions presented in Sec. 4.2.6. As discussed in Sec. 4.2.6 the additional gluon contributions induce several problems. The first problem can be clearly seen in both plots in Fig. 12. The Goldstone nature of the DM candidate χ requires that the SI cross section scales  with the corresponding DM mass m χ [7], implying that the SI cross section vanishes in the zero DM mass limit, since the Goldstone nature of the DM candidate is restored. Note that this particular behaviour is only expected for vanishing momentum transfer as assumed in the calculation. Our approach (neglecting the gluon contributions) shows for both benchmark points (left and right in Fig. 12) the desired behaviour for small DM masses m χ which does not happen when the additional gluon contributions to the SI cross section are included. As for the second problem related to the approximation performed in the two-loop diagrams, it is not clear how it would reflect on the results. What we can see from the plots is that for large DM masses both approaches yield similar results. The difference is roughly a factor three induced by the inclusion of the gluon contributions. Further, the results presented in Ref. [9] are exactly reproduced only if we include the additional gluon contributions. The important point here is to understand that unless a complete 2-loop calculation of the gluon contribution is performed, nothing can be said about the inclusion of approximate calculations of some diagrams.
We calculated all contributing diagrams in the general R ξ gauge in order to be able to check for missing gauge cancellations. As it turned out, our result is completely gauge independent. For the proper cancellation of all gauge dependencies the Goldstone triangle diagrams in Fig. 7 were crucial. They were needed to properly cancel the gauge dependencies introduced in the vertex corrections. These diagrams are often overseen in the literature. However, the inclusion of the additional gluon contributions introduces gauge dependencies which are not cancelled. We define the relative change where σ indicates the SI cross section calculated with the additional gluon contributions in the general R ξ gauge and in the Feynman gauge (ξ = 1), respectively. In the left plot of Fig. 13 we show the results for the relative change as a function of the gauge parameter ξ. The color scheme follows that of Fig. 12. Again the benchmark scenario presented in Tab. 1 is used to determine the SI cross section. Obviously, when the additional gluon contributions are included, the variation of the gauge parameter ξ changes the SI cross section significantly preventing to make reasonable predictions for the NLO SI cross section. Hence, not only the correct DM mass dependence is lost with the inclusion of the additional gluon diagrams but also a strong gauge dependence is introduced.
Finally we will discuss the contribution of the third diagram in Fig. 8. In the right plot of Fig. 13 we show the relative difference where σ| nogb is the SI cross section calculated without the effective two-loop vertex (third diagram in Fig. 8) and σ the SI cross section as presented. We varied the DM mass m χ while keeping the other input parameters (same as in Tab. 1) fixed. The most dominant effect of the gluon boxes are obtained for small or large DM masses. Despite that, the overall impact given by the gluon boxes is in the sub-percentage region. Hence, not taking into account the gluon box diagrams and thereby treating all diagrams with external gluons in Fig. 8 consistently would not significantly alter the overall result. In addition to the phenomenological discussion of the SI cross section of the PNGDM we implemented the model in BSMPT [57,58] allowing us to check for a first-order electroweak phase transition in the early universe. For simplicity we force the vacuum expectation value of the DM field component χ in Eq. (2) to be equal to zero at all temperatures while determining the global NLO minimum at finite temperature. This way we ensure the stability of the DM candidate and its DM nature. All parameter points above the neutrino floor provide an NLO stable vacuum in the sense that the vacuum ground state of the one-loop effective potential (at zero temperature) is the same as the tree-level ground state. However, all parameter points provide a weak phase transition v c /T c < 1, where v c is the SU (2) VEV (v) at the critical temperature T c . The critical temperature is defined as the temperature where the one-loop effective potential has two degenerate minima. A more involved study in the phase structure at finite temperature of the PNGDM might enable a strong first order electroweak phase transition. For instance, allowing the DM field component to evolve a non-zero VEV at finite temperature leads to interesting phenomenological consequences. These studies are left for future work.

Conclusions
In this work we have calculated the NLO corrections to the spin independent scattering cross section of a scalar DM particle off a nucleon in a Pseudo Nambu-Goldstone DM model. This model has a scalar potential invariant under a global U (1) symmetry softly broken such that a pseudo Nambu-Goldstone boson originates from the broken symmetry. The cross section was first shown to be proportional to the Dark Matter velocity in Ref. [6]. Therefore there was the need to perform the calculation at NLO. There were two independent calculations that appeared very close in time [7,9].
The first calculation [7] was performed by considering from the effective Lagrangian only the first term q C q S O q S . Instead of nuclear matrix elements for the proton an effective Higgs-nucleon coupling was used. Because in this case the one loop result for the Wilson coefficient is independent of the quark masses, it factorises, and it turns out that the Higgs-nucleon effective coupling is the sum of the nuclear matrix elements. This calculation reproduces the correct dependence of the cross section in the limit of vanishing Dark Matter and is at least 90% of the total cross section depending on the parameter points. Hence, relative to this work we have now included the terms C g S O g S + q C q T O q T . In the second calculation that appeared in Ref. [9] all terms in the effective Lagrangian were used. As previously discussed, the problem in this calculation resides in the gluon diagram contributions (specifically the first two diagrams in Fig. 8). The 2-loop diagram is not effectively calculated and instead an approximation is performed such that a proper matching between the heavy quarks and the gluon operators cannot be performed. The only way to solve the problem would be to perform the complete 2-loop calculation. The approximation leads therefore to the fact that the Goldstone nature of the DM candidate is not recovered and that the these gluon contributions present a strong gauge dependence contrary to the rest of the calculation. Therefore these additional gluon contributions should be dropped unless the full two-loop calculation is performed.
It is also worth mentioning that although we have used a different renormalisation scheme than the ones from the two previous calculations our results show a very similar behaviour when the different contributions are compared. Finally, we showed that with the present constraints most of the allowed points are below the neutrino floor and only experiments in the far future will be able to probe them.