Electroweak Corrections to Dark Matter Direct Detection in a Vector Dark Matter Model

Although many astrophysical and cosmological observations point towards the existence of Dark Matter (DM), the nature of the DM particle has not been clarified to date. In this paper, we investigate a minimal model with a vector DM (VDM) candidate. Within this model, we compute the cross section for the scattering of the VDM particle with a nucleon. We provide the next-to-leading order (NLO) cross section for the direct detection of the DM particle. Subsequently, we study the phenomenological implications of the NLO corrections, in particular with respect to the sensitivity of the direct detection DM experiments. We further investigate more theoretical questions such as the gauge dependence of the results and the remaining theoretical uncertainties due to the applied approximations.


Introduction
The first indirect hints of the existence of Dark Matter (DM) were reported more than 80 years ago [1] (see [2] for a historical account). Over the years, confirmations from different sources have established the existence of DM. Still, as of today, these are all indirect evidences and all based on gravitational effects. Therefore, it may come as a surprise that today the discussion on the properties of DM is about particles which is a consequence of the findings from astronomy and cosmology. Direct, indirect and collider searches for DM refer in most cases to particles belonging to some extension of the Standard Model (SM), and all experimental data from the different sources favour a weakly interacting massive particle (WIMP) with a velocity of the order of 200 km/s. That is, DM is non-relativistic.
Many experiments have been proposed for the direct detection of DM on Earth. It was shown in Ref. [3] that DM particles that undergo coherent scattering with nuclei, i.e. spin-independent scattering, are the ones with larger scattering rates, and therefore they can be detected more easily. The scattering rates depend on the material of the detector, on the underlying cosmological model through the assumption of an approximately constant DM density, ρ 0 = 0.3 GeV/cm 3 , on the DM velocity and finally on the DM-nucleon cross section. Although there are uncertainties associated with the determination of all the parameters, the need for an increased precision in the DM-nucleon cross section calculation has led several groups to invest in the calculation of higher-order corrections, both strong and electroweak, to the scattering cross section [4][5][6][7][8][9][10][11][12][13].
Although the hypothesis of DM as a particle is now the strongest and most intensely studied conjecture to explain the data, there are no hints on the exact nature of the particle itself. Among the several possibilities, in this study we will focus on a minimal model with a vector DM candidate. The model is a very simple extension of the SM where a dark vector χ µ with a gauged U (1) χ symmetry and a complex SM-gauge singlet S are added to the SM field content. We end up with a vector DM (VDM) candidate χ µ and a new CP-even scalar that mixes with the SM scalar field coming from the doublet.
The electroweak corrections to the coherent scattering of the DM candidate χ µ first require the renormalisation of the VDM model and second, the extraction of the spin-independent contributions from the loop corrections to the effective couplings of the Lagrangian, L eff , which couple two DM particles and two quarks. These will then constitute the corrections to the tree-level effective couplings from L eff .
The paper is organised as follows: in Section 2 we present the VDM model and in Section 3 we describe its renormalisation. In Section 4 we discuss the scattering of scalar DM off nuclei at leading order (LO) while in Section 5 we calculate the electroweak corrections to the cross section. In Section 6 we present and discuss our results. In the conclusions, Section 7, we summarise our findings. Feynman rules and technical details are left to the appendices.

The Vector Dark Matter Model
The VDM model discussed in this work is an extension of the SM, where a complex SM-gauge singlet S is added to the SM field content [14][15][16][17][18][19][20]. The model has a new U (1) χ gauged symmetry, under which solely the gauge singlet S is charged. As the symmetry is gauged, a new vector boson appears in the theory, which is denoted by χ µ .
In order to obtain a stable VDM candidate we assume a Z 2 symmetry. The dark gauge boson χ µ and the scalar field S transform under the Z 2 symmetry as follows χ µ → −χ µ and S → S * , (2.1) and the SM particles are all even under Z 2 , which precludes kinetic mixing between the gauge bosons from U (1) χ and the SM U (1) Y . As the singlet S is charged under the dark U (1) χ , its covariant derivative reads D µ S = (∂ µ + i g χ χ µ ) S , (2.2) where g χ is the gauge coupling of the dark gauge boson χ µ .
The most general Higgs potential invariant under the SM and the Z 2 symmetries can be written as in terms of the squared mass parameters µ 2 H , µ 2 S and the quartic couplings λ H , λ S and κ. The neutral component of the Higgs doublet H and the real part of the singlet field each acquire a vacuum expectation value (VEV) v and v S , respectively. The expansions around their VEVs can be written as where Φ H and Φ S denote the CP-even field components of H and S. The CP-odd field components σ H and σ S do not acquire VEVs and are therefore identified with the neutral SM-like Goldstone boson G 0 and the Goldstone boson G χ for the gauge boson χ µ , respectively, while G ± are the Goldstone bosons of the W bosons. The minimum conditions of the potential yield the tadpole equations which allow the scalar mass matrix to be expressed as (2.7) The treatment of the tadpole contributions in the mass matrix will be discussed in Section 3 while describing the renormalisation of the tadpoles. The mass eigenstates h 1 and h 2 are obtained through the rotation with the orthogonal matrix R α as The diagonalisation of the mass matrix yields the mass values m h 1 and m h 2 of the two scalar mass eigenstates. The mass of the VDM particle will be denoted as m χ . The parameters of the potential Eq. (2.3) can then be expressed in terms of the physical parameters using the relations The SM VEV v ≈ 246 GeV is fixed by the W boson mass. The mixing angle α can be chosen without loss of generality to be α ∈ − π 2 , π 2 . (2.14) The requirement of the potential to be bounded from below is translated into the following conditions

Renormalisation of the VDM Model
In order to calculate the electroweak (EW) corrections to the scattering process of the VDM particle with a nucleon we need to renormalise the VDM model. There are four new independent parameters relative to the SM that need to be renormalised. We choose them to be the non-SM-like scalar mass, m h 2 , the rotation angle α, the coupling g χ and the DM mass m χ . 1 In the following, we will present the renormalisation of the VDM model including the gauge and Higgs sectors.
Having chosen the complete set of free parameters in the theory, we start by replacing the bare parameters p 0 with the renormalised ones p according to (3.16) where δp is the counterterm for the parameter p. Denoting a general scalar or vector field by Ψ, the renormalised field is expressed in terms of the field renormalisation constant Z Ψ as where Ψ 0 stands for the bare and Ψ for the renormalised field, respectively. In case of mixing field components, √ Z Ψ is a matrix.
Gauge Sector: Since we have an extended gauge sector compared to the SM we will give all counterterms explicitly. Due to the imposed Z 2 symmetry under which only the dark gauge boson χ µ is odd, kinetic mixing between the gauge bosons of the U (1) χ and to U (1) Y is not possible. This means that there is no interaction between the gauge sector of the SM and the new dark gauge sector. Since this symmetry is only broken spontaneously, gauge bosons from the two sectors will not mix at any order of perturbation theory and therefore the field renormalisation constants are defined independently in each sector. We choose to renormalise the theory in the mass basis. The replacement of the parameters in the two gauge sectors reads e → e + δZ e e , (3.18d) where m W and m Z are the masses of the electroweak charged and neutral gauge bosons W ± and Z, respectively, e is the electric coupling constant, and g the weak SU (2) coupling. The gauge boson fields are renormalised by their field renormalisation constants δZ, The on-shell (OS) conditions yield the following expressions for the mass counterterms of the gauge sector 20) where T denotes the transverse part of the self-energies Σ ii (i = W, Z, χ). Expressing the electric charge in terms of the Weinberg angle θ W e = g sin θ W , with and using OS conditions for the renormalisation of the electric charge allows for the determination of the counterterm δg in terms of the mass counterterms δm W , δm Z and δZ e 2 δZ e = 1 2 The wave function renormalisation constants guaranteeing the correct OS properties are given by
As for the gauge coupling from the dark sector, g χ , since there is no obvious physical quantity to fix the renormalisation constant, we will renormalise it using the MS scheme, which will be described in detail in Section 3.1.
Higgs Sector: In the VDM model we have two scalar fields which mix, namely the real component Φ H of the Higgs doublet and the real component Φ S of the singlet, yielding the mass eigenstates h 1 and h 2 . This mixing has to be accounted for in the field renormalisation constants (see Eq. (3.17)) so that the corresponding matrix reads In the mass eigenbasis, the mass matrix in Eq. (2.7) yields The tadpole terms in the tree-level mass matrix are bare parameters. At next-to-leading order (NLO) they obtain a shift that corresponds to the change of the vacuum state of the potential through electroweak corrections. To avoid such vacuum shifts at NLO, we renormalise the tadpoles such that the VEV remains at its tree-level value also at NLO. This requires the introduction of tadpole counterterms δT i such that the one-loop renormalised one-pointT i function vanishesT Since we formulate all counterterms in the mass basis it is convenient to rotate the tadpole parameters in their corresponding mass basis as well, using the same rotation matrix R α , For the mass counterterms of the Higgs sector we replace the mass matrix as with the one-loop counterterm (3.31) In Eq. (3.31) we neglect all terms of order O (δαδT i ) since they are formally of two-loop order. Using OS conditions and Eq. (3.31) finally yields the following relations for the counterterms (i = 1, 2) 3.1 Renormalisation of the Dark Gauge Coupling g χ As previously mentioned, the dark gauge coupling g χ cannot be linked to a physical observable, which prevents the usage of OS conditions for its renormalisation. Therefore the coupling will be renormalised using the MS scheme. As the UV divergence is universal, we just need a vertex involving g χ . We choose the triple h 1 h 1 h 1 vertex. First we write where A VC stands for the amplitude for the virtual corrections to the vertex and A CT is the amplitude for the vertex counterterm. We will henceforth drop the index h 1 h 1 h 1 for better readability. The counterterm amplitude A CT consists of two contributions, and The trilinear Higgs self-coupling reads (expressing v through 2m W /g) The divergent part of δg χ is then given by In Fig. 1 we present the set of diagrams used to calculate A VC . The one-loop diagrams were generated with FeynArts [21] for which the model file was obtained with SARAH [22][23][24][25] and the program package FeynCalc [26,27] was used to reduce the amplitudes to Passarino-Veltmann integrals [28]. The numerical evaluation of the integrals was done by Collier [29][30][31][32]. The counterterm g χ in the MS scheme is then obtained as where γ E denotes the Euler-Mascheroni constant.

Renormalisation of the Scalar Mixing Angle α
The final parameter that needs to be renormalised is the mixing angle α. Again, this is a quantity that cannot be related directly to an observable, except if we would use a process-dependent renormalisation scheme which is known to lead to unphysically large counterterms [33]. The renormalisation of the mixing angles in SM extensions was thoroughly discussed in [33][34][35][36][37][38][39][40][41][42][43][44]. In this work we will use the KOSY scheme, proposed in [45,46], which connects for the derivation of the angle counterterm the usual OS conditions of the scalar field with the relations between the gauge basis and the mass basis. The bare parameter expressed through the renormalised one and the counterterm reads 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 With Eq. (3.34) the mixing angle counterterm reads In our numerical analysis we will use two more renormalisation schemes for δα: the MS scheme and a process-dependent scheme. In the MS scheme we only take the counterterm δα into account in the divergent parts in D = 4 dimensions. Applying dimensional regularisation [47,48], these are the terms proportional to 1/ , where D = 4 − 2 . Both the KOSY scheme and the MS scheme lead to a gauge-parameter dependent definition of δα This is avoided if δα is defined through a physical process.
In our process-dependent renormalisation scheme for α, discussed in the numerical results, we define the counterterm δα through the process h → τ τ , where h denotes the SM-like scalar of the two h i (i = 1, 2). The counterterm is defined by requiring the NLO decay width to be equal to the LO one. The NLO corrections involve infrared (IR) divergences stemming from the QED corrections. Since they form a UV-finite subset, this allows us to apply the renormalisation condition solely on the weak sector thus avoiding the IR divergences, i.e. we require for the NLO and LO amplitudes of the decay process where 'weak' refers to the weak part of the NLO amplitude. The h coupling to ττ depends on the mixing angle α as

51) and the LO amplitude reads
with u(p τ ) (ū(p τ )) denoting the spinor (anti-spinor) of the τ with four-momentum p τ . Dividing the weak NLO amplitude into the LO amplitude, the weak virtual corrections to the amplitude, and the corresponding counterterm part, and we get the mixing angle counterterm in the process-dependent scheme as Here A ct δα=0 denotes the complete counterterm amplitude but without the contribution from δα.

Dark Matter Direct Detection at Tree Level
In the following we want to set our notation and conventions used in the calculation of the spin-independent (SI) cross section of DM-nucleon scattering. The interaction between the DM and the nucleon is described in terms of effective coupling constants. The major contribution to the cross section comes from the light quarks q = u, d, s and gluons. For the VDM model the effective operator basis contributing to the SI cross section is given by [49] where G a µν (a = 1, ..., 8) denotes the gluon field strength tensor and O q µν the quark twist-2 operator corresponding to the traceless part of the energy-momentum tensor of the nucleon [50,51], Operators suppressed by the DM velocities and the momentum transfer of the DM particle to the nucleon are neglected in the analysis. Furthermore, we neglect contributions introduced by the gluon twist-2 operator O g µν , since these contributions are one order higher in the strong coupling constant α s [49].
For vanishing momentum transfer and on-shell nucleon states, the nucleon matrix elements are given by where N denotes a nucleon, N = p, n, and m N is the nucleon mass. Furthermore, q N (2),q N (2) are the second moments of the parton distribution functions of the quark q(x) and the antiquark q(x), respectively. The four-momentum of the nucleon is denoted by p. The numerical values for the matrix elements f N Tq , f N T G and the second moments q N (2) andq N (2) are given in App. A. The SI effective coupling of the VDM particle with the nucleons is obtained from the nucleon expectation value of the effective Lagrangian, Eq. (4.56), by applying Eqs. (4.59), which yields In the contribution from the quark twist-2 operator all quarks below the energy scale ∼ 1 GeV have to be included, i.e. all quarks but the top quark. The SI scattering cross section between the VDM particle and a nucleon, proton or neutron (N = p, n), is given by Note that the sum in the first term of Eq. (4.60) only extends over the light quarks. The leading-order gluon interaction with two VDM particles is mediated by one of the two Higgs bosons which couple to the external gluons through a heavy quark triangle diagram, cf. Fig. 2. The charm, bottom and top quark masses are larger than the energy scale relevant for DM direct detection and should therefore be integrated out for the description of the interaction at the level of the nucleon. By calculating the heavy quark triangle diagrams and then integrating out the heavy quarks we obtain the related operator in the heavy quark limit. This is equivalent to calculating the diagram in Fig. 3 with heavy quarks Q = c, b, t, and replacing the resulting tensor structure m QQ Q with the effective gluon operator as follows [12,13,52] corresponding to the effective leading-order VDM-gluon interaction in Eq. (4.57).
For the tree-level contribution to the SI cross section the t-channel diagrams depicted in Fig. 3 have to be calculated for vanishing momentum transfer. The respective Wilson coefficient for the effective operator in Eq. (4.56) is extracted by projecting onto the corresponding tensor structure, m q qq. Accounting for the additional symmetry factor of the amplitude, this yields finally the following f q factor for the quarks, As explained above, the heavy quarks Q = b, c, t contribute to the effective gluon interaction. By using Eq. (4.62), the Wilson coefficient for the gluon interaction, f G , can be expressed in terms of f q for q = c, b, t, resulting in the SI LO cross section The twist-2 operator does not contribute at LO. The obtained result is in agreement with Ref. [20] 3 .

Dark Matter Direct Detection at One-Loop Order
As a next step, we want to include the NLO EW corrections in the calculation of the SI cross section. For this, we evaluate the one-loop contributions to the Wilson coefficients f q and f G in front of the operators in Eq. (4.57). At this order, also the Wilson coefficient g q is non-zero. The additional topologies contributing at NLO EW are depicted in Fig. 4. Note that we do not include vertex corrections to the h iq q vertex. They are partly given by the nuclear matrix elements and beyond the scope of our study. For the purpose of our investigation, we assume them to be encoded in the effective coupling factors of the respective nuclear matrix elements.
In the following, we present the calculation of each topology separately.

Vertex Corrections χχh i
The effective one-loop coupling χχh i is extracted by considering loop corrections to the coupling χχh i , where we take the DM particles to be on-shell and assume a vanishing momentum for the Higgs boson h i . The amplitude for the NLO vertex including the polarisation vectors ε ( * ) of the external VDM particles, is given by with the leading-order amplitude i A LO χχh i , the virtual vertex corrections i A VC χχh i and the vertex counterterm i A CT χχh i . Denoting by p the four-momentum of the incoming VDM particle, the tree-level amplitude is given by The vertex counterterm amplitudes for i = 1, 2 read with the counterterms δg χχh i (i = 1, 2) for the couplings g χχh 1 = 2g χ m χ sin α (5.69) Assuming zero momentum transfer is equivalent to assuming p in = p out . Note that this assumption is stricter than simply assuming p 2 in = p 2 out , since this only implies the same masses for the incoming and outgoing particles. The additional new tensor structure (denoted by ∼ NLO) is given by (5.72) The additional NLO tensor structure vanishes by assuming p in = p out , and because for freely propagating gauge bosons we have ε(p) · p = 0. The counterterms in Eq. (5.68) cancel all UV-poles of the virtual vertex corrections in Fig. 5 which has been checked both analytically and numerically. Accounting for the symmetry factor of the amplitude and projecting onto the corresponding tensor structure, the vertex corrections are plugged in the generic diagram in Fig. 4(a) which contributes to the operator χ µ χ µ m qq q. We will refer to the resulting contribution as f vertex q . Since the expression it quite lengthy, we do not give the explicit formula here.

Mediator Corrections
We proceed in a similar way for the mediator corrections. We calculate the self-energy corrections to the two-point functions with all possible combinations of external Higgs fields and plug these into the one-loop propagator in the generic amplitude in Fig. 4(b). The self-energy contribution to the h i h j propagator (i, j = 1, 2) reads with the renormalised self-energy matrix where the mass matrix M and the tadpole counterterm matrix δT are defined in Eq. (3.27). The Z-factor matrix δZ corresponds to the matrix with the components δZ h i h j defined in Eq. (3.34). Projecting the resulting one-loop correction on the corresponding tensor structure, we obtain the effective one-loop correction to the Wilson coefficient of the operator χ µ χ µ m qq q induced by the mediator corrections as with the rotation matrix R α defined in Eq. (2.8).

Box Corrections
We now turn to the box corrections. The generic set of diagrams representative of the box topology is depicted in Fig. 6. In the following, we present the treatment of box diagrams contributing to the SI cross section. In order to extract for the spin-independent cross section the relevant tensor structures from the box diagram, we expand the loop diagrams in terms of the momenta p q of the external quark that is not relativistic [12]. Since we are considering zero momentum transfer, the incoming and outgoing momenta of the quark are the same, Note that as in the case of the vertex corrections this requirement is stricter than requiring that the squared momenta are the same, since this only implies same masses for incoming and outgoing particles. Assuming small quark momenta, and because the mass of the light quarks is much smaller than the energy scale of the interaction, allows for the simplification of the propagator terms arising in the box diagrams through the expansion, where l is the loop momentum of the box diagram, m q the mass of the quark and where we use m 2 q = p 2 q . After applying this expansion to the box diagrams, the result has to be projected onto the required tensor structures contributing to the operators in Eq. (4.57). The box diagrams contribute to X µ X µ m qq q and the twist-2 operators. By rewriting [13,50,51] Figure 7: The full two-loop gluon interaction with the DM candidate (left) and the effective two-loop interaction after integration out the heavy quarks (right). and therefore can be dropped. We refer to these one-loop contributions to the corresponding tree-level Wilson coefficients as f box q and g box q . As discussed in Refs. [12,13] the box diagrams also induce additional contributions to the effective gluon interaction with the VDM particle that have to be taken into account in the Wilson coefficient f G in Eq. (4.57b). The naive approach of using the same replacement as in Eq. (4.62) to obtain the gluon interaction induces large errors [12]. To circumvent the overestimation of the gluon interaction without performing the full two-loop calculation, we adopt the ansatz proposed in Ref. [13]. For heavy quarks compared to the mediator mass, it is possible to derive an effective coupling between two Higgs bosons and the gluon fields. Using the Fock-Schwinger gauge allows us to express the gluon fields in terms of the field strength tensor G a µν , simplifying the extraction of the effective two-loop contribution to f G . Integrating out the top-quark yields the following effective two-Higgs-two-gluon coupling [13] where the effective coupling d eff G of Ref. [13] has to be adopted to our model. First of all we only have scalar-type mediators, given by the Higgs bosons h i , so that the mixing angle φ SM of Ref. [13] which quantifies the CP-odd admixture, is set to Second, the coupling of the Higgs bosons h i to the top quark differs depending on which Higgs boson is coupled, so that the effective coupling in Eq. (5.79) becomes with the rotation matrix R α defined in Eq. (2.8). The effective coupling allows for the calculation of the box-type diagram in Fig. 7 (right).
In Ref. [13], the full two-loop calculation was performed. The comparison with the complete two-loop result showed very good agreement between the approximate and the exact result for mediator masses below m t . Our model is structurally not different in the sense that the mediator coupling to the DM particle (a fermion in Ref. [13]) is also a scalar particle so that the results obtained in Ref. [13] should be applicable to our model as well. Moreover, the box contribution to the NLO SI direct detection cross section is only minor as we verified explicitly.
The diagram in Fig. 7 (right) yields the following contribution to the Lagrangian where C ij denotes the contribution from the triangle loop built up by h i , h j and the VDM particle. It has to be extracted from the calculated amplitude of Fig. 7 (right). Using Eq. (4.57b) the contributions by the box topology to the gluon-DM interaction are given by

The SI One-Loop Cross Section
In the last sections we discussed the extraction of the one-loop effective form factors for the operators in Eq. (4.57). The NLO EW SI cross section can then be obtained by using the effective one-loop form factor with the Wilson coefficients at one-loop level given by Like at LO, the heavy quark contributions of f vertex q and f med q have to be attributed to the effective gluon interaction. With the LO form factor given by where f LO q has been given in Eq. (4.63), we have for the NLO EW SI cross section at leading order in α S , (5.87)

Numerical Analysis
In our numerical analysis we use parameter points that are compatible with current theoretical and experimental constraints. These are obtained by performing a scan in the parameter space of the model and by checking each data set for compatibility with the constraints. In order to do so, the VDM model was implemented in the code ScannerS [53, 54] which automatises the parameter scan. We require the SM-like Higgs boson (denoted by h in the following) to have a mass of m h = 125.09 GeV [55]. With ScannerS, we check if the minimum of the potential is the global one and if the generated points satisfy the theoretical constraints of boundedness from below and perturbative unitarity. We furthermore impose the perturbativity constraint g 2 χ < 4π. Furthermore, the model has to comply with the experimental Higgs data. In the VDM model, the Higgs couplings to the SM particles are modified by a common factor given in terms of the mixing angle α, that is hence constrained by the combined values for the signal strengths [55]. Through an interface with HiggsBounds [56][57][58] we additionally check for collider bounds from LEP, Tevatron and the LHC. We require agreement with the exclusion limits derived for the non-SM-like Higgs boson at 95% confidence level. Among these searches the most stringent bound arises from the search for heavy ZZ resonances [59]. Still, the bounds for the mixing angle α derived from the measurement of the Higgs couplings are by far the most relevant. In order to check for the constraints from the Higgs data, the Higgs decay widths and branching ratios were calculated with sHDECAY [54] 5 , which includes the state-of-the-art higher-order QCD corrections. The code sHDECAY is based on the implementation of the models in HDECAY [60,61].
Concerning the DM constraints, information on the DM particle from LHC searches through the invisible width of the SM Higgs boson were taken into account [56][57][58]. Furthermore, the DM relic abundance has been calculated with MicrOMEGAs [62][63][64][65], and compared with the current experimental result from the Planck Collaboration [66], We do not force the DM relic abundance to be in this interval, but rather require the calculated abundance to be equal to or smaller than the observed one. Hence, we allow the DM to not saturate the relic density and therefore define a DM fraction where (Ωh 2 ) χ stands for the calculated DM relic abundance of the VDM model. DM indirect detection also provides constraints on the VDM model. The annihilation into visible states, mainly into ZZ, W + W − , bb and light quark pairs, can be measured by Planck [66], if it manifests itself in anisotropies of the cosmic microwave background (CMB), by Fermi-LAT [67] if it comes form the γ-ray signals in the spheroidal dwarf galaxies, and by AMS-02 [68,69] if it originates from e ± excesses in the Milky Way. As shown in Ref. [70], the Fermi-LAT upper bound on the DM annihilations is the most stringent one. In order to obtain the bound, we follow Ref. [67] in their claim that all final states give approximately the same upper bound on the DM annihilation cross sections. Hence we use the Fermi-LAT bound from Ref. [67] on bb when m χ m b , and on light quarks for m χ < m b . In the comparison with the data, the DM fraction in Eq. (6.89) has to be taken into account, and an effective DM annihilation cross section is defined by with f χχ and σ χχ computed by MicrOMEGAs.
The sample was generated taking into account the experimental bounds on the DM nucleon SI cross section at LO. The most stringent bound on this cross section is the one from the XENON1T [71,72]   The ranges of the input parameters of the scan performed to generate viable parameter sets are listed in Table 1. From here on, we denote the non-SM-like of the two CP-even Higgs bosons m h i (i = 1, 2) by m φ , the SM-like Higgs boson is called m h . Note, that in ScannerS the scan is performed over m χ and v S instead of m χ and g χ . The corresponding g χ values are given by g χ = m χ /v S . Only points with g 2 χ ≤ 4π are retained. Note that we vary α in the range [−π/4, π/4] to optimize the scan. This is possible due to the bound on sin α that comes from the combined signal strength measurements of the production and decay of the SM-like Higgs boson [55]. The remaining input parameters, gauge, lepton and quark masses, electric coupling, Weinberg angle and weak SU (2)

Results
In the following we present the LO and NLO results for the spin-independent direct detection cross section of the VDM model. We investigate the size of the NLO corrections and their phenomenological impact. We furthermore discuss the gauge dependence of the results and the influence of the renormalisation scheme on the NLO results. If not stated otherwise, results are presented for the Feynman gauge, i.e. the gauge parameter ξ 6 is set equal to one, ξ = 1. In the NLO results, the default renormalisation scheme for the mixing angle α is the KOSY scheme, cf. Subsection 3.2.

The SI Direct Detection Cross Section at Leading Order
In Fig. 8 we show in grey the LO results of the direct detection cross section for all points of the VDM model that are compatible with our applied constraints, as a function of the DM mass m χ . Note, that we also include the perturbativity limit on g χ , g 2 χ < 4π. The result is compared to the Xenon limit shown in blue. Note that, in order to be able to compare with the Xenon limit, we applied the correction factor f χχ to the LO and NLO direct detection cross section, cf. Eq. (6.91). Since the compatibility with the Xenon limit is already included in the selection of valid parameter points, all cross section values lie below the blue line (modulo the size of the grey points). As can be inferred from Fig. 8, the LO cross section can be substantially smaller than the present sensitivity of the Xenon experiment, by more than 10 orders of magnitude.

Results for m φ < m t
We now investigate the dependence of the LO and NLO direct detection cross section on g χ and the size of the NLO corrections for the parameter sets featuring a non-SM-like Higgs boson with a mass m φ < m t . For these, the approximate treatment of the NLO box contributions discussed in Subsection 5.3 can be applied. In Fig. 9 we display for all parameter sets passing our constraints that additionally feature m φ < m t the LO direct detection cross section in the left panel and the NLO result in the right panel, as a function of m φ . The color code quantifies the coupling g χ ≤ √ 4π. Note, that here and in the following we do not apply the correction factor f χχ Eq. (6.89) on the direct detection limit, as long as we do not directly compare with the Xenon limit. This is why the LO cross section in Fig. 9 is larger than in Fig. 8.
Both the LO and the NLO contribution to the SI direct detection cross section are proportional to f LO q and therefore proportional to g χ , sin 2α and (m 2 . This behaviour is reflected in Fig. 9. We observe that the LO cross section increases with g χ , more specifically g 2 χ (yellow points) and drops for m φ = m h = 125.09 GeV. The NLO corrections on the other hand increase with g 3 χ . The reason is that, as we explicitly verified, the NLO corrections are dominated by the vertex corrections. The vertex corrections are proportional to g 2 χ , so that the NLO contribution to the cross section scales as 2 Re(f LO q f vertex * q ) ∝ g 3 χ , in contrast to the LO contribution that is proportional to g 2 χ . In total the K-factor, i.e. the ratio between NLO and LO cross section, therefore increases with g χ .
Being proportional to f LO q the NLO corrected cross section also drops for m φ = m h , so that the sensitivity of the direct detection experiment is not increased after inclusion of the NLO corrections; the blind spots remain also at NLO. In our scan we furthermore did not find any parameter points where a specific parameter combination leads to an accidental suppression at LO that is removed at NLO. There is a further blind spot when α = 0. However, in this case the SM-like Higgs boson has exactly SM-like couplings and the new scalar decouples from all SM particles except for the coupling with the SM-like Higgs boson. In this scenario the SM-like Higgs decouples from the vector dark matter particle, and, depending on the mass of the second scalar and its coupling strength with the SM-like Higgs boson, we may end up with two dark matter candidates with the second scalar being metastable. The study of such a scenario is beyond the scope of this paper and we do not consider this case further here. The K-factor is depicted in Fig. 10, as a function of m φ (left) and σ LO (right). The colour code indicates the size of g χ . The K-factor is mostly positive and the bulk of K-facture values ranges between 1 and about 2.3. As mentioned above, the K-factor increases with g χ , as can also be inferred from the figure, in particular from Fig. 10 (right).
In this and all other plots, we excluded points with m φ ≈ m h and K-factors where |K| > 2.5. We found that for m φ ≈ m h the interference effects between the h and φ contributions, that become relevant here, largely increase the (dominant) vertex contribution f vertex q to the effective NLO form factor. It exceeds by far the LO form factor f LO q . Depending on the sign of f vertex q , the NLO cross section, which is proportional to 2 Re(f LO q f vertex * q ), is largely increased or suppressed, inducing for large negative NLO amplitudes negative NLO cross sections. In these regions, the NLO results are therefore no longer reliable. Two-loop contributions might lead to a better perturbative convergence, but are beyond the scope of this paper. We deliberately did not take into account one-loop squared terms to remove the negative cross sections. Such an approach would only include parts of the two-loop contributions. Whether or not they approximate the total two-loop results well enough can only be judged after performing the complete two-loop calculation. This is why we chose the conservative approach to exclude these points from our analysis.
In Fig. 11, we show the K-factor as function of σ LO , but with the colour code indicating the size of sin 2 2α (left) and m χ right. There is no clear correlation between the K-factor and sin 2 2α or m χ . These plots furthermore show, that there is no correlation between the maximum size of σ LO and m χ or sin 2 2α, while the maximum σ LO values require large g χ values, cf. Fig. 10 (right).

Results for m φ > m t
We now turn to the parameter region of our sample of valid points where the approximation described in Subsection 5.3 is a priori not valid. We cannot judge the goodness of the approximation in this parameter region without doing the full two-loop calculation which is beyond the scope of this paper. We can check, however, if we see some unusual behaviour compared to the results for parameter sets with m φ < m t , where the approximation can be applied.  Figure 12 shows the K-factor as a function of the LO SI direct detection cross section. The size of g χ is indicated by the color code. We only take into account parameter samples compatible with all constraints and where m φ > m t . As already observed and discussed for the parameter sample with m φ < m t , also here the K-factor increases with g χ . Overall, the bulk of points reaches larger K-factors than for m φ < m t but remains below 2.5. So, the behaviour of the K-factor does not substantially differ from the results for m φ < m t . The comparison of the approximate and exact result in Ref. [13] showed that the difference in the box contribution between the two results does not exceed one order of magnitude for a pseudoscalar mediator with mass 1 TeV 7 and remains small even for scalar mediator masses up to 1 TeV (cf. Fig. 4 in Ref. [13]). Together with the fact that the box contribution makes up only for a small part of the NLO SI direct detection cross section 8 , we conclude that our approximate NLO results for parameter sets with larger mediator masses should also be applicable in this parameter region.

Gauge Dependence
As has been discussed in Ref. [33] for the 2HDM and in Ref. [37] for the N2HDM the renormalisation of the mixing angle α in the KOSY scheme leads to gauge parameter dependent results. We therefore check here the gauge dependence of our NLO results by performing the calculation in the general R ξ gauge and comparing it with our default result in the Feynman gauge ξ = 1.
We introduce the relative gauge dependence ∆ ξ σ NLO , defined as where σ NLO ξ denotes the NLO SI direct detection cross section calculated in the general R ξ gauge and σ NLO ξ=1 the result in the Feynman gauge. In Fig. 13 we show ∆ ξ σ NLO as a function of the gauge parameter ξ, for two sample parameter points of our valid parameter set, called 7 We estimate this by extrapolating Fig. 4 (left) in Ref. [13] to 1 TeV 8 We explicitly verified that the box form factor f box q remains below the vertex correction form factor f vertex q . In particular, for K-factors above 1, the box form factor remains more than two orders of magnitude below the vertex form factor.  As can be inferred from Fig. 13, we clearly see a gauge dependence of the NLO results. The relative gauge dependence is, however, small with values below the few per cent level for a rather large range of ξ variation. Note also, that a gauge parameter dependence a priori is no problem as long as it is made sure that the explicit value of the gauge parameter is accounted for when interpreting the results.

Renormalisation Scheme Dependence
We now investigate the influence of the renormalisation scheme and scale on the NLO result. For this, we show in Fig. 14 the K-factor for the whole data sample passing our constraints for three different renormalisation schemes of the mixing angle α as a function of the LO cross section. The chosen schemes have been described in Subsection 3.2 and are the KOSY scheme (yellow points), the process-dependent scheme (green) and the MS scheme (violet). The scale applied in the MS scheme is µ 0 = 1 GeV. The KOSY scheme has been shown to lead to a gaugeparameter dependent definition of the counterterm δα [33,37]. This is also the case for the MS scheme. As can be inferred from the plot, the MS scheme additionally leads to unnaturally large NLO corrections with K-factor values up to about 10 8 for our data sample (not shown in the plot). This has been known already from previous investigations in the 2HDM [33] and N2HDM [37]. The process-dependent scheme has the virtue of implying a manifestly gaugeindependent definition of the mixing angle counterterm. However, also here the NLO corrections are unacceptably large with values up to about 10 9 , so that also this scheme turns out to be unsuitable for practical use. This behaviour has also been observed in our previous works [33,37].
We therefore conclude that the KOSY scheme should be used in the computation of the NLO corrections. The fact that it is gauge dependent is no problem as long as the chosen gauge is clearly stated when presenting results. Moreover, by applying a pinched scheme, the gauge dependence can be avoided, cf. Ref. [33]. This is left for future work.  Figure 14: The K-factor versus the LO SI direct detection cross section for the whole data sample passing all constraints and for three different renormalisation schemes of α: the KOSY scheme (yellow), the process-dependent scheme (green), the MS scheme (violet).
The uncertainty due to missing higher-order corrections can be estimated by varying the renormalisation scheme or by varying the renormalisation scale. The comparison of the KOSY with the other two renormalisation schemes makes no sense as the latter lead to unacceptably large corrections. The KOSY scheme does not allow us to vary the renormalisation scale, so that we cannot provide an estimate of the uncertainty due to missing higher order corrections. We conclude with the remark that the variation of the renormalisation scale between 1/2 and 2 times the scale µ 0 in the MS scheme leads to a variation of the NLO cross section of about 16% -in contrast to the unphysically large corrections that are to be traced back to the blowing-up of the MS counterterm of α.

Phenomenological Impact of the NLO Corrections on the Xenon Limit
We now turn to the discussion of the phenomenological impact of our NLO results. In Fig. 15 (left) we show the LO direct detection cross section (blue points) and the NLO result (orange) compared to the Xenon limit (blue-dashed), as a function of the DM particle mass. For the consistent comparison with the Xenon limit we applied the correction factor f χχ (Eq. 6.91) to the LO and NLO cross section in both plots of Fig. 15. In the left figure we plot all parameter points where the LO cross section does not exceed the Xenon limit but the NLO result does. This plot shows that for a sizeable number of parameter points, the compatibility with the experimental constraints would not hold at NLO any more. This demonstrates that the NLO corrections are important and need to be accounted for in order to make reliable predictions about the viable parameter space of the VDM model.
In the right plot we display the same quantities, but only for parameter points of our data  This implies we only consider parameter points where the LO cross section is much smaller than the Xenon limit, but the NLO cross section is of the order of the Xenon limit. We learn from this figure that although LO results might suggest that the Xenon experiment is not sensitive to the model, this statement does not hold any more when NLO corrections are taken into account. These results confirm the importance of the NLO corrections when interpreting the data.

Conclusions
In this paper, we investigated a minimal model with a VDM particle. We computed the NLO corrections to the direct detection cross section for the scattering of the VDM particle off a nucleon. We developed the renormalisation of the model, proposing several renormalisation schemes for the mixing angle α of the two physical scalars of the model. We computed the leading corrections, including relevant two-loop box contributions to the effective gluon interaction in the heavy quark approximation. With the box contributions to the NLO cross section being two orders of magnitude below the leading vertex corrections, we estimated the error induced by the approximation to be small. Interference effects of the two scalar particles that become important for degenerate mass values on the other hand, were found to be large and require further investigations beyond the scope of this paper, namely the computation of the complete two-loop contributions. Outside this region, the perturbative series is well-behaved and Kfactors of up to about 2.5 were found.
We further investigated the impact of the chosen renormalisation scheme for α. While the process-dependent renormalisation of α is manifestly gauge-parameter independent, it was found to lead to unphysically large corrections. This did not improve by choosing the gaugeparameter dependent MS scheme. A renormalisation scheme exploiting the OS conditions of the scalar fields on the other hand, leads to moderate K-factors, while being manifestly gaugeparameter dependent. For the proper interpretation of the data, therefore, the choice of the gauge parameter has to be specified here.
We found that the NLO corrections can either enhance or suppress the cross section. With K-factors of up to about 2.5, they are important for the correct interpretation of the viability of the VDM model based on the experimental limits on the direct detection cross section. The NLO corrections can increase the LO results to values where the Xenon experiment becomes sensitive to the model, or to values where the model is even excluded due to cross sections above the Xenon limit. In case of suppression, parameter points that might be rejected at LO may render the model viable when NLO corrections are included.
The next steps would be to investigate in greater detail the interesting region of degenerate scalar masses and study its implication on phenomenology in order to further be able to delineate the viability of this simple SM extension in providing a VDM candidate.