Loop-induced direct detection signatures from CP-violating scalar mediators

We investigate direct detection signatures of dark matter particles interacting with quarks via a light spin-0 mediator with general CP phases. Since tree-level scattering may be strongly suppressed in the non-relativistic limit, loop contributions play an important role and can lead to observable signals in near-future experiments. We study the phenomenology of different mediator masses and CP phases with an emphasis on scenarios with maximal CP violation and Higgs portal models. Intriguingly, the sum of the rates obtained at tree- and loop-level can give a characteristic recoil spectrum not obtainable from a single type of interaction. We furthermore develop a novel method for decomposing the two-loop contribution to effective interactions between dark matter and gluons into two separate one-loop diagrams, which in our case substantially simplifies the calculation of the important top-quark contribution.


Introduction
Experiments aiming to directly detect the interactions of dark matter (DM) particles in underground laboratories have made tremendous progress over the past decades and place some of the strongest bounds on the parameter space of many DM models [1]. Indeed, these experiments have become so sensitive that they can be relevant even for DM models where the leading order interactions are momentum-or velocity-suppressed [2]. As a result there has been a rapidly growing interest in the general effective field theory (EFT) of non-relativistic interactions between DM and nuclei [3][4][5][6][7][8][9][10][11][12]. In these models it becomes essential to include loop effects, which may reintroduce spin-independent interactions and thereby substantially boost the expected event rates [13][14][15][16][17][18][19].
Particular attention has been paid to models in which DM scattering is mediated by a pseudoscalar exchange particle [19][20][21][22][23], motivated partially by the interesting implications for collider [24][25][26][27][28][29][30][31][32][33] and flavour [34][35][36][37][38][39] physics. At leading order the resulting interactions are so strongly suppressed in the non-relativistic limit that they are well below the -1 -JHEP06(2019)052 "neutrino floor" which indicates the ultimate reach of direct detection experiments [40]. However, several recent studies have shown that loop-induced spin-independent interactions can change this picture dramatically, in particular when taking into account the interactions between the pseudoscalar mediator and the SM Higgs boson required by gauge invariance [41][42][43][44]. In fact, ref. [44] pointed out that for this particular model even twoloop processes give a relevant contribution and need to be properly included for an accurate estimate of experimental sensitivities.
In the present work we generalise these results by considering spin-0 mediators that couple to DM and Standard Model (SM) quarks with arbitrary CP phases. We furthermore treat the coupling between the mediator and SM Higgs bosons as a free parameter and thus remain agnostic about the underlying ultraviolet (UV) completion. A particular emphasis is placed on the impact of two-loop processes. We show that, at least for heavy quarks, accurate results can be obtained by first integrating out the heavy quark and then performing all further calculations in the resulting EFT. This approach substantially simplifies and speeds up the evaluation of direct detection constraints.
We find that for general CP phases loop-induced spin-independent interactions may be strong enough to lead to detectable signals in near-future direct detection experiments, such as LZ [45] or XENONnT [46]. The importance of our results are illustrated for a number of relevant scenarios. We show that for DM models with maximal CP violation (as studied e.g. in the context of self-interacting DM [47]) loop effects can be comparable to the leadingorder contribution and change the shape of the recoil spectrum in important ways. Large effects are also found in the CP-violating Higgs portal model, which has been the subject of several recent studies [48][49][50]. In both cases loop-induced interactions enable direct detection experiments to probe parameter regions that would otherwise be out of reach.
The paper is structured as follows. In section 2 we briefly introduce the general model with free CP phases and then present our central results on how to perform the mapping onto the low-energy EFT relevant for DM direct detection. We discuss in detail the importance of two-loop processes and the matching onto non-relativistic effective operators. Specific applications of the general results are presented in section 3, where we also calculate the sensitivity of present and future direct detection experiments. We summarise our findings and conclude in section 4. Detailed results from our one-loop and two-loop calculations are presented in the appendices A and B, respectively. Finally, appendix C provides details on nuclear form factors.

Loop effects in direct detection
We investigate a simplified model of a Dirac fermion DM particle χ interacting with SM fermions f through a general spin-0 mediator a with mass m a greater than the bottomquark mass m b : where φ χ and φ SM are CP phases, v ≈ 246 GeV is the electroweak vacuum expectation value, m f are the SM fermion masses and g χ as well as g SM denote the couplings of a to -2 -

JHEP06(2019)052
DM and SM fermions, respectively. We have further assumed Yukawa-like couplings in agreement with the hypothesis of minimal flavour violation (MFV) [51] such that flavour physics constraints on the universal coupling g SM are weakened (see section 3.2). 1 For φ χ = φ SM = 0 we recover the well-known simplified model of a scalar mediator, whereas for φ χ = φ SM = π/2 we obtain a CP-conserving theory with a pseudoscalar mediator [52].
In the former case constraints on the model from direct detection experiments are very strong, whereas in the latter case they are almost entirely absent [41,42,44]. Here we will treat the CP phases as free parameters in order to study the impact of different phase combinations on the predictions for direct detection experiments. The simplified model in eq. (2.1) does not respect all gauge symmetries of the SM before electroweak symmetry breaking. The interactions between a and SM fermions are therefore expected not to appear in isolation but in combination with additional interactions between a and the SM Higgs boson h. In the present work, we will not discuss how these different interactions can be linked in specific UV completions. Instead, we introduce an additional free parameter λ ah and supplement eq. (2.1) by the interaction term We will show that this interaction can play a relevant role in the phenomenology of this model. Moreover, it will be of particular importance in section 3.3 where we will identify a with the SM Higgs boson h itself. Note that we neglect additional interaction terms involving two Higgs bosons. Although such terms are in general expected to be present, they do not give any relevant contribution to the calculation of direct detection signatures. We finally note that for φ SM = π/2 the mediator a can mix with the SM Higgs boson, giving rise to direct interactions of the SM Higgs boson with DM particles. This mixing is however required to be small given that the observed Higgs behaves SM-like in current experiments. Furthermore, the Higgs boson mass is much larger than the values of m a that we will consider, such that its contribution to direct detection is suppressed [42,44]. We will therefore not consider Higgs mixing within this work but emphasize that it would be straightforward to include these contributions using the results presented below.

Low-energy effective Lagrangian
To calculate event rates in direct detection experiments, we need to determine the effective interactions between DM and quarks that result from the three types of diagrams illustrated in figure 1. For the discussion below it will be useful to distinguish between interactions that lead to spin-independent (SI) and to spin-dependent (SD) scattering in the non-relativistic limit. 2 Starting with the tree-level exchange of a illustrated in the left panel of figure 1, 1 In a generic MFV scenario a slightly more general Lagrangian than eq. (2.1) can be written down, as different couplings to up-and down-type quarks are allowed. For the scope of this work, however, we will focus on the case of one universal coupling. 2 Note that here and below we use the term "spin-independent" to refer to all types of interactions that do not depend on the nucleus spin, irrespective of whether or not they are suppressed in the non-relativistic limit. Accordingly, the term "spin-dependent" refers to all interactions that are not spin-independent, including momentum-suppressed interactions. Indeed, unsuppressed spin-dependent interactions are absent in the model that we consider. we obtain where the sum runs over all quark species. Here we have defined the tree-level coefficient and have kept the dependence on the two CP phases explicit. Next we consider the Higgs-mediated exchange shown in the middle panel of figure 1, which maps onto the purely spin-independent interaction where the sum again runs over all quarks. We have further introduced the triangle coefficients 8) in terms of the loop functions C 0 (m 2 χ , m 2 a , m 2 χ ) and C 2 (m 2 χ , m 2 a , m 2 χ ), which are given in appendix A.1.
Finally, we have to take into account the box diagram in the right panel of figure 1. We expand the amplitude in terms of the quark momentum, which is the smallest scale in the diagram [44], and obtain Computational details and the expressions of the different box diagram coefficients C box i,q are given in appendix A.2. Note that all of these coefficients share a common factor of g 2 χ g 2 SM m 2 q /v 2 , which also constitutes the only quark dependence. In eq. (2.9) we have also introduced the twist-2 quark operator Since the corresponding form factors are evaluated at the scale of the Z boson mass m Z , we include the charm and bottom quark in the corresponding sums in eq. (2.9) [54,55]. However, none of the heavy quarks have been included in the remaining terms of eqs. (2.9) and (2.10), because they require a different treatment, which will be discussed next.

Effective description of two-loop processes
As the charm, bottom and top quark are heavier than the energy scale relevant for DM direct detection experiments, they should be integrated out of the theory aiming to describe interactions at the level of nuclei. For the tree-level and Higgs-induced triangle diagram this can be done simply by replacing the heavy quarks by the corresponding effective gluon interaction obtained from triangular heavy-quark loops [56] where G aµν is the gluon field strength tensor and G aµν = 1 2 µναβ G a αβ with the convention 0123 = 1. This procedure is justified for these two diagrams since the two steps of integrating out the mediator a and integrating out the heavy quarks factorise.
The situation is however very different for the box diagram in the right panel of figure 1. In this case one cannot make a simple factorization argument to integrate out heavy quarks. This is visualised in figure 2, which shows the two-loop diagrams that need to be computed to obtain the effective interactions between DM and gluons. Any attempt to simplify this calculation by first integrating out the mediator a and then using eq. (2.12) would neglect the contribution from the diagram on the right. For m Q m a , m χ the twoloop computation hence cannot be simplified in this way without introducing potentially  Figure 3. Illustration of the decomposition of the two-loop process for m Q m a , m χ . After first integrating out the heavy quark Q (first arrow) one can then match the resulting one-loop diagram onto effective DM-gluon interactions (second arrow). The black dot represents an effective interaction corresponding to a higher-dimensional operator.
large errors [44]. In the opposite case of m Q m a , m χ it was argued in ref.
[44] that a simplification is not possible because one cannot expand the box diagram amplitude in terms of the external quark momentum, which is no longer the smallest scale in the diagram. It was in particular stressed that for the top quark a full two-loop computation is mandatory.
As we are now going to demonstrate, however, for m Q m a , m χ it is in fact possible to decompose the underlying two-loop process into two separate one-loop diagrams by integrating out the heavy quark Q first and the mediator a afterwards. This approach, in which no diagrams are neglected, is illustrated in figure 3. Provided the mediator is light compared to the heavy quark, it is thus possible to simplify the calculations significantly.
In the following we will be mostly interested in the case m a m t , such that the approach outlined above can be applied to the top quark. Therefore, we first consider the loop involving the top quark separately and integrate it out by performing a 1/m t expansion of the (in total six) corresponding amplitudes. We employ Package-X [57] for the evaluation and expansion of the loop computations. This then yields the following leading order effective Lagrangian coupling a to gluons (2.14) Here we have included a symmetry factor of 1/2 and defined 3 which are both independent of the top-quark mass. Now performing the second step 3 Note that d eff G vanishes for certain values of φSM such that one would need to include higher orders. However, these specific cases are not of interest in the present work. While d eff G also vanishes for specific values of φSM, the same is true for the full expression d full G , see eq. (B.31) in appendix B.2, i.e. this is not a result of the heavy quark expansion.  visualised in figure 3, we obtain for the effective two-loop approach where the effective two-loop coefficients read An analogous calculation for the bottom and charm quark only gives a useful approximation if m a m c , m b . For heavier mediator masses it is in general unavoidable to perform the full two-loop calculation to accurately estimate the corresponding contributions (see appendix B.2 for more details). However, for the specific coupling structure that we are interested in, bottom and charm quark are found to give only a small contribution. 4 It is hence possible to obtain a very good approximate result of the total heavy quark contribution to the effective DM-gluon interactions by including only the top-quark contribution using our effective approach. This is illustrated in figure 4, where we plot the absolute value of the coefficient C G,S as a function of the mediator mass (left panel) and of the DM mass m χ (right panel). The effective approach for the top quark (indicated by the solid red line) and the corresponding two-loop calculation (dotted green) show very good agreement for m a m t across the whole range of DM masses. Including also bottom and charm quark in the two-loop JHEP06(2019)052 calculation has only slight influences for small values of m a (dashed grey). Similar results can be obtained for the other coefficients. 5 We conclude that it is possible to simplify the full two-loop calculation in the case of m Q m a , m χ , which makes it possible to circumvent the full two-loop calculation entirely if the top quark is expected to give the dominant contribution. We will therefore use the effective approach for the remainder of this work.

Matching onto effective operators
In this section we match the effective interactions of DM with quarks and gluons onto non-relativistic DM-nucleon interactions in order to obtain predictions for direct detection experiments. The first step is to perform the matching of quark and gluon fields onto nucleon fields, which yields the following effective Lagrangian: (2.20) where N = p, n is a nucleon field and 'CPV' indicates terms that only arise when CP is violated. The coefficients C eff depend on the various coefficients we derived in the previous two sections as well as on the nuclear form factors that parametrise the quark and gluon contents of a nucleon. Note that in general the nuclear form factors and hence the effective coefficients are different for protons and neutrons: C eff,p = C eff,n . For the SI coefficients, we find as well as where the nuclear form factors f N q,G , q N (2) andq N (2) are defined in appendix C. 5 For specific parameter points cancellations might occur within the coefficients CG and C G like in d eff G for φSM ≈ π/4. In this parameter region the two-loop result and the effective approach differ. However, this discrepancy does not affect any of the scenarios studied in detail below.

JHEP06(2019)052
Likewise, we obtain for the SD coefficients The form factors F q/N P and F N G are given in appendix C. Because of non-negligible contributions from the π and η pole, these form factors depend on the momentum exchange q µ between DM and nucleons.
In the non-relativistic limit the effective Lagrangian from eq. (2.20) can be matched onto a basis of effective operators: where the operators O N i depend only on the spins S χ and S N of DM and the nucleon, respectively, as well as on the momentum transfer q and the DM-nucleon relative velocity v [3,4,58]. For the model that we consider, only four different operators are generated, namely (2.26) The corresponding coefficients can be directly read off from L eff χN [4]: Note that like the form factors F q/N P and F N G the coefficients c N 6 and c N 10 also depend on the momentum transfer. This final step completes the derivation of the effective interactions relevant for DM direct detection from the general Lagrangian of a spin-0 mediator given in eq. (2.1).

JHEP06(2019)052 3 Phenomenological implication
In this section we use the results from above to predict the differential event rates in past and future direct detection experiments and to calculate the resulting exclusion limits and expected sensitivities. In models that predict dominantly spin-independent scattering, this can be done by simply calculating the corresponding scattering cross section is the DM-nucleon reduced mass. For c p 1 ≈ c n 1 the differential event rate with respect to recoil energy E R is then simply given by where ρ 0 is the local DM density, A is the mass number of the target nucleus and F 2 (E R ) denotes the nuclear form factor. The factor g(v min ) = v min f (v)/v dv denotes the velocity integral as a function of the minimum velocity v min (E R ) = m A E R /(2 µ 2 ) with m A being the mass of the target nucleus and µ being the corresponding reduced mass. Direct detection experiments typically assume this particular form of the differential cross section in order to produce exclusion limits and quote expected sensitivities in terms of σ SI p as a function of m χ .
In the presence of additional interactions, however, the calculation of the differential event rate becomes much more involved. We do not review the corresponding formalism here and instead refer to refs. [3,4,10]. Crucially, for momentum-dependent interactions it is no longer possible to capture the model prediction in terms of a single cross section at fixed momentum transfer which can then be compared to published exclusion limits. To evaluate experimental sensitivity it thus becomes necessary to reproduce experimental analyses for the appropriate recoil spectra and include information on detection efficiencies and background levels in order to obtain approximate likelihood functions. This process has been automated for the most general set of non-relativistic effective operators in the public code DDCalc v2.1 [49,59], which includes an extensive database of existing and planned direct detection experiments. Furthermore, DDCalc contains an automated interface with DirectDM [60], which we use for the matching of the spin-dependent coefficients in eq. (2.27) and the evaluation of the corresponding nuclear form factors. We can therefore simply pass the coefficients C eff,N calculated for our model to DDCalc and obtain the likelihoods for existing direct detection experiments and the predicted number of events in future experiments. In the following, we will indicate the regions of parameter space that are excluded by the most recent XENON1T results [61] and the regions that predict at least 5 events in the next-generation LZ experiment [45]. 6 Similar exclusion limits are obtained from the Panda-X [2,62] and LUX [63,64] experiments, while comparable sensitivities are expected for the XENONnT experiment [46].  Figure 5. Direct detection constraints as a function of the CP-violating phases φ χ and φ SM . For both panels we have fixed g χ = g SM = 1 and λ ah = 0. The dotted lines indicate the ratio of the predicted number of events in LZ when including both tree-level and loop-level diagrams and when only including tree-level diagrams. This ratio can be smaller than unity due to destructive interference.

General CP phases
We first visualize the general aspects of our model by considering in figure 5 the most general case, in which φ χ and φ SM can take arbitrary values between 0 (corresponding to purely scalar couplings) and π/2 (corresponding to purely pseudoscalar couplings). For the purpose of this figure we have fixed g SM = g χ = 1 and λ ah = 0 and consider two different combinations of m χ and m a in the two panels. Note that we assume that the correct relic density is reproduced at each point in both plots without invoking a specific mechanism. The blue shading indicates the parameter region excluded by XENON1T, while the dashed green line provides an estimate for the reach of LZ. The black dotted lines indicate the ratio of the total number of predicted events in LZ to the number of events predicted at tree-level.
For φ χ , φ SM π/2 the tree-level exchange of a dominates the spin-independent coefficient C SI eff,N from eq. (2.20) and therefore also the whole scattering process. In such a scenario current direct detection bounds rule out a large part of the parameter space and constrain g SM to be very small [47,65]. As the two phases approach π/2, tree-level scattering becomes more and more suppressed, leading to a reduced sensitivity of direct detection experiments and a greater importance of loop effects.
For φ χ = 0 and φ SM = π/2, i.e. the top-left corner of figure 5, CP violation is maximal. In this case the tree-level contribution maps onto the non-relativistic operator O N 10 , which is suppressed in the non-relativistic limit and furthermore depends on the spin of the nucleus. Existing direct detection constraints can thus be evaded even with O(1) couplings [21]. However, spin-independent contributions arise at loop-level and can dominate the event rate and yield potentially observable signals. The importance of loop-effects can also be seen for φ χ ≈ 0 and 10 −4 π/2 − φ SM 10 −3 , where the total event rate is smaller than the one predicted at tree-level due to the destructive interference between spin-independent -11 -JHEP06(2019)052 interactions present at tree-level and those induced at loop-level. We discuss the case of maximal CP violation in more detail in section 3.2.
For the opposite scenario of φ χ = π/2 and φ SM = 0, i.e. the bottom-right corner of figure 5, the tree-level contribution to spin-independent scattering maps onto the nonrelativistic operator O N 11 , which depends on the DM spin and the momentum transfer. While the scattering cross section does receive a coherent enhancement in this case, it is suppressed by an additional factor of m 2 N /m 2 χ . We will therefore study the influence of purely spin-independent contributions emerging at loop-level in the context of the CPviolating Higgs-portal model in section 3.3.
Finally, in the top-right corner of figure 5, corresponding to almost purely pseudoscalar interactions, the loop-induced event rate dominates over the tree-level prediction by many orders of magnitude. However, as observed previously [44], the sensitivity of direct detection experiments is strongly suppressed in this limit, so that the case of pure pseudoscalar interactions is out of reach for current direct detection experiments. A crucial conclusion from figure 5 is that loop effects become increasingly important as experimental sensitivity improves. For the couplings and masses considered, XENON1T is only sensitive to those regions in parameter space where loop-induced interactions give a sub-leading contribution. LZ on the other hand will be sensitive to interactions that are more strongly suppressed at tree-level, giving greater importance to an accurate calculation of loop-level contributions.

Maximal CP violation
Let us take a closer look at the case φ χ = 0 and φ SM = π/2, corresponding to the topleft corner in figure 5. In this case spin-independent interactions are completely absent at tree-level, making loop effects particularly important. Indeed, for the masses and couplings considered in figure 5 this scenario is not excluded by the bounds from XENON1T but can be tested with LZ. However, the loop contributions depend sensitively on the strength of the couplings, which enter quadratically into the Wilson coefficients. In order to fully assess the importance of loop effects, it is therefore important to consider alternative constraints on the couplings g SM , g χ and λ ah .
For given values of m a , m χ and g SM we can fix g χ by the requirement that the observed DM relic abundance can be explained in terms of thermal freeze-out via the annihilation processes χχ → qq and χχ → aa. If the latter process is kinematically allowed, i.e. for m a < m χ , it will typically give the dominant contribution for g SM 1, such that the required value for g χ becomes independent of g SM . In this limit, we find g χ ∝ m 1/2 χ with g χ = 1 for m χ ≈ 500 GeV. For larger g SM the calculation becomes more involved and we use micrOmegas v5.0.6 [66] to determine the required value for g χ numerically.
The coupling of the light spin-0 boson to SM particles can be constrained through a range of flavour physics observables. For m a m B ≈ 5.2 GeV, these constraints are very strong and effectively exclude the possibility of obtaining observable direct detection signatures [37]. However, almost all of these constraints disappear for larger values of m a . Bounds from radiative Υ decays [67,68] extend to slightly larger masses, but also disappear for m a 7 GeV. Provided the pseudoscalar couples also to leptons (with coupling strength g SM m /v), another important constraint arises from B s → µ + µ − , which can arise -12 -

JHEP06(2019)052
from loop-induced flavour-changing interactions with an off-shell mediator. The resulting branching ratio is given by [37,69] where C SM 10 = −4.103 and Λ is the scale of new physics (such as additional charged Higgs bosons needed in a gauge-invariant UV completion). For Λ = 1 TeV and assuming m a m Bs , this expression simplifies to The branching ratio of B s → µ + µ − has been measured with a precision of 20% [70] and is found to be in agreement with the SM prediction [71]. To obtain an approximate bound on g SM we therefore require the new-physics contribution not to exceed 40% of the SM value. This gives g SM 1.0 m a 10 GeV . (3.5) In other words, even for spin-0 bosons as light as 10 GeV the coupling strength g SM can be of order unity. Constraints of comparable strength have been derived from LHCb dark photon searches within a di-muon channel, see refs. [32,72]. 7 The situation is quite different for the coupling λ ah between a and the SM Higgs boson. This coupling induces the decay h → aa with partial width [48] The presence of this decay mode gives rise to exotic Higgs decays and leads to a suppression of the Higgs signal strength in the conventional channels. While the former provide a promising strategy for future searches [32], at present the strongest constraints come from a global fit of the measured properties of the SM-like Higgs boson at ATLAS and CMS [78]. These fits imply BR(h → aa) < 0.34, corresponding to Γ h→aa 2 MeV, when simultaneously allowing for modifications of the Higgs boson production cross section, or BR(h → aa) < 0.13, corresponding to Γ h→aa 0.6 MeV, when assuming the production cross section to be given by the SM prediction. For m a m h /2, these bounds translate to λ ah 0.02 and λ ah 0.01, respectively. We will conservatively show the weaker bound in the following.  Figure 6. Constraints on g SM (left) and λ ah (right) as a function of m χ in a model with maximal CP violation (φ χ = 0, φ SM = π/2). At each point the coupling g χ is fixed in such a way that the observed DM relic abundance is reproduced. The dotted lines indicate the ratio of the predicted number of events in LZ from loop-induced spin-independent interactions and from tree-level momentumsuppressed interactions. green line. Dotted black lines in the left panel indicate the ratio of loop-induced spinindependent interactions and tree-level momentum suppressed interactions in terms of the number of predicted events in LZ. As expected, the importance of loop effects grows with increasing g SM and with increasing m χ , corresponding to increasing g χ . The kinks for m χ ≈ 175 GeV result from the fact that for larger DM masses annihilation into top quarks becomes kinematically allowed and provides an efficient annihilation channel, reducing the required value of g χ .
Since direct constraints on g SM are quite weak, we find large regions of parameter space where the model can be discovered by LZ. If the interactions of DM arise dominantly from λ ah , on the other hand, the strong constraints from Higgs measurements imply that there remains only a small region of allowed parameter space that can be explored with LZ. We note that the constraints in the right panel are completely independent of φ SM and would hence also apply to a pure pseudoscalar.
For parameter points close to the XENON1T exclusion bound in the left panel loop effects give a sizeable contribution to the total event rate in direct detection experiments. This observation is illustrated further in figure 7, which compares the predicted differential event rates at tree-level and loop-level in LZ for m χ = 200 GeV, m a = 15 GeV and g SM = 0.6, corresponding to g χ = 0.6. The tree-level interactions are momentum-suppressed and therefore vanish in the limit E R → 0, leading to a maximum around E R ∼ 30 keV. The differential event rate from loop-induced spin-independent interactions, on the other hand, decreases monotonically with increasing recoil energy. Intriguingly, the two contributions conspire to give a total event rate that is approximately constant across the entire search region. Such a spectrum cannot be obtained from any single non-relativistic operator and could therefore, given enough statistics, be used to identify models like the one discussed here. A similar interplay between tree level and loop level can arise for φ χ = π/2, φ SM = 0, in which case the tree-level process is coherently enhanced but suppressed by a factor m N /m χ in c 11 , see eq. (2.27). The two scenarios however differ in their dependence on the target material. In particular, if tree-level scattering is spin-dependent, it will be absent in target materials with no nuclear spin, leading to a monotonically falling recoil spectrum from loop-induced spin-independent interactions.
Let us finally revisit the discussion of how to approximate two-loop effects in our model. We compare in figure 8 the spin-independent scattering cross section obtained with our approach (outlined in section 2.2) with the result of a full two-loop calculation including all heavy quarks. The left panel corresponds to the case of maximal CP violation (φ χ = 0, φ SM = π/2), the right panel corresponds to the pure pseudoscalar case (φ χ = φ SM = π/2). In both cases we fix g χ by the relic density requirement and set g SM = 0.6, consistent with the bounds discussed above (which are independent of φ χ ). We find very good agreement between the two approaches, confirming our approach for integrating out top quarks and neglecting the contribution from bottom and charm quarks. In the right panel we also show the cross section obtained if the pseudoscalar is integrated out before all heavy quarks, as previously suggested in refs. [41,42]. 8 As pointed out previously [44], this approach leads to a vast overestimation of the loop contribution.

CP-violating Higgs portal
As a final example for the importance of loop-effects we consider the fermionic Higgs portal model [48][49][50]:  In the left panel we consider maximal CP violation with φ SM = π/2 and φ χ = 0, whereas in the right panel we fix φ SM = φ χ = π/2, i.e. pure pseudoscalar phases, and also show the curve corresponding to previous calculational approaches of the two-loop diagram. In both panels λ ah is set to zero and g χ is fixed such that the correct relic density is reproduced. Note that additional contributions to the differential event rate from momentum-dependent interactions at tree-level may lead to stronger exclusion limits than the ones shown in this plot. discussed so far. After electroweak symmetry breaking, however, the following interactions are generated: where λ denotes the quartic Higgs self-coupling and (3.10) We can therefore directly apply all the results from section 2 with the replacements The factor of 6 in the last expression is necessary to ensure that the correct Feynman rule is obtained in spite of different combinatorial factors. The free parameters of this model are hence m χ , λ hχ /Λ and φ. Note that, since we study loop processes within this model, one should in principle include all operators involving DM particles and SM fields that contribute at the order of 1/Λ 2 , in particular dimension six operators coupling the DM vector and axial-vector current to the corresponding SM quark currents. Here we implicitly assume the absence of new spin-1 particles at the high energy scale Λ that would induce such operators. Operators including scalar, pseudoscalar or tensor couplings between DM and quarks would generally be accompanied with a factor of m q and would therefore only contribute at higher order, i.e. 1/Λ 3 .
For φ = 0 the model violates CP and spin-independent scattering is suppressed proportional to cos 2 φ. As φ approaches π/2, loop effects are therefore expected to become increasingly important. We confirm this expectation in figure 9, which shows constraints on λ hχ /Λ as a function of m χ . Dotted lines indicate the ratio of loop-induced spin-independent interactions to tree-level momentum-suppressed interactions (in terms of the expected number of events in LZ). In the parameter range that can be probed by direct detection experiments, this ratio is significantly larger than unity, implying that the sensitivity of direct detection experiments stems almost exclusively from loop-induced interactions. 9 In figure 9 we also indicate the parameter regions excluded by the constraint BR(h → inv) < 0.26 [79,80] as well as the combinations of λ hχ /Λ and m χ for which the observed DM relic abundance can be reproduced via annihilations into SM particles [48]. The requirement of EFT validity, λ hχ /Λ < 2π/m χ [49], is satisfied in the entire parameter region shown in figure 9. We find that for φ = π/2 constraints from direct detection experiments are rather weak and only probe parameter regions where the standard freeze-out calculation predicts χ to be a sub-dominant DM component. For these parameter regions we

JHEP06(2019)052
implicitly assume that the abundance of χ is set by a non-standard mechanism (e.g. a particle-antiparticle asymmetry) such that χ accounts for all of the DM. If, on the other hand, bounds from direct detection experiments are rescaled based on the abundance of χ obtained from standard freeze-out, as done e.g. in ref. [49], loop-induced direct detection signals do not provide relevant constraints on the CP-violating Higgs portal model for the foreseeable future.
We point out that within this model two additional loop diagrams with one insertion each of hχiγ 5 χ and h 2χ iγ 5 χ contribute to the amplitudes relevant for direct detection. These diagrams are however UV divergent, which indicates a dependence on the specific UV completion of the effective Higgs portal operator. Replacing the UV divergence with log (Λ 2 /m 2 χ ) and setting Λ = 1 TeV, we find that these loops can be numerically important and increase the predicted event rates (see appendix D for additional details). Nevertheless, these additional contributions are still not large enough for near-future direct detection experiments to reach the relic density line shown in figure 9. To make this statement more precise would require the choice of a specific UV completion.

Conclusions
Future direct detection experiments will reach such a high level of sensitivity to the interactions between DM and quarks that loop effects become increasingly important. This is particularly true in models where tree-level scattering is suppressed, such that loop-induced interactions may give the dominant contribution and yield potentially observable signals. In the present work we have studied such a set-up in the context of a spin-0 particle a mediating the interaction between DM and SM fermions. In contrast to previous studies, we allow general CP phases and therefore cover scalar, pseudoscalar and CP-violating interactions. Moreover, we include a trilinear coupling between a and the SM Higgs boson which generally arises in UV completions of this model and can have important phenomenological consequences.
For certain combinations of CP phases standard spin-independent contributions are strongly suppressed or even fully absent at tree-level, such that a proper calculation of the interactions induced at loop-level is crucial. In our model, these arise from Higgs-induced triangle diagrams, box diagrams for light quarks (both shown in figure 1) as well as the twoloop process involving heavy quarks shown in figure 3. In particular the two-loop process gives an important contribution, which is difficult to estimate without performing the full calculation. To address this challenge, we have presented a novel approach for simplifying the two-loop calculation significantly for heavy quark masses (schematically illustrated in figure 3). Provided the top quark gives the dominant contribution and the mediator is light compared to the top quark, this approach makes it possible to circumvent the two-loop calculation entirely and obtain an accurate estimate that is much easier to calculate and implement. A comparison between the two approaches is provided in figure 4.
As illustrated in figure 5, loop effects are most important when at least one of the CP phases is close to π/2 (corresponding to pseudoscalar interactions). Moreover, they gain in importance as the sensitivity of direct detection experiments improves. A particu--18 -

JHEP06(2019)052
larly interesting observation is that the recoil rates induced at tree-and loop-level can be comparable, resulting in a roughly constant event rate over the whole energy window (see figure 7). Since such a spectrum cannot be generated from a single type of interaction, it will be very interesting to perform a detailed statistical analysis of how to discriminate the model studied here from alternative hypotheses.
Finally, we have studied the impact of spin-independent loop-induced interactions on the CP-violating fermionic Higgs portal model. Our results show that loop-level effects allow future direct detection experiments to probe parameter regions that would be otherwise inaccessible. Nevertheless, loop-level contributions are still too small to enable direct detection experiments to reach the parameter regions preferred by thermal freeze-out if the CP phase is close to π/2.
Based on the results presented in this work, we conclude that a general spin-0 mediator offers an interesting possibility to evade current direct detection bounds even with O(1) couplings while still maintaining promising detection prospects for future years. It will therefore be important to investigate how such a simplified model can arise from a more complete theory, such as an extended Higgs sector with spontaneous CP breaking. Such an embedding will provide new insights on the relations between the different couplings and allow for a more accurate analysis of the constraints from flavour physics and precision observables.

Acknowledgments
We thank Giorgio Arcadi and Sebastian Wild for discussions and Joachim Brod for valuable comments on the manuscript. This work is funded by the Deutsche Forschungsgemeinschaft (DFG) through the Emmy Noether Grant No. KA 4662/1-1 and the Collaborative Research Center TRR 257 "Particle Physics Phenomenology after the Higgs Discovery".

A One-loop Wilson coefficients
In this appendix we provide details on the one-loop calculations relevant for section 2.1.

A.1 Loop functions
We define the Passarino-Veltman functions C i that appear in our one-loop calculations according to the standard notation [81] The X 2 and Y 2 functions are given by [82]

JHEP06(2019)052
which will reappear in the full two-loop approach in appendix B.1. Finally, we define the Z functions All of these functions can be calculated readily with Package-X [57]. As examples, we quote the expression for C 0 (p 2 , m 2 , M 2 ), and Z 11 (M 2 , M 2 , m 2 ) evaluated using the on-shell condition p 2 = M 2 , The remaining coefficients can be computed analogously.

A.2 Box diagram computation and coefficients
For the computation of the box and its crossed diagram shown in figure 10 we follow the procedure from ref. [44], which allows us to derive the coefficients for the twist-2 operators. We first start with the amplitude, which can be expressed as where the crossed diagram is obtained by the replacement k → −k withinū q (p q ) . . . u q (p q ) and we have suppressed the sum over the quark species. Now we expand the amplitude in terms of p q , as this is the smallest scale involved in the diagram: where we have employed the on-shell condition (p q ) 2 = m 2 q . The amplitude then reads We can now identify the loop functions defined in appendix A.1, construct the corresponding effective Lagrangian and use the following decomposition: This then yields the effective box diagram Lagrangian given in eqs. (2.9) and (2.10) with the coefficients given by where we have used the shorthand notation

B Details on two-loop calculations
In this appendix we provide details on the two-loop calculations relevant for section 2.2.

B.1 Loop functions
For the two-loop computation presented in appendix B.2, we will need further loop functions. The Passarino-Veltman functions D i read in their standard notation [81] We further define the X n and Y n functions by [82] Using partial fraction decomposition, we can derive the following relations for X n where the loop functions denoted by B 0 (p 2 , m 2 , M 2 ) and B 1 (p 2 , m 2 , M 2 ) are implemented in LoopTools [83]. 10 We further need the derivatives of these functions with respect to m 2 . (B.16) 10 We could in principle evaluate the functions Xn and Yn directly with Package-X [57] and then perform the numerical integration appearing in the two-loop calculation using their explicit expressions. Numerical stability improves, however, when the functions are decomposed as presented here. The two-loop computation was recently presented for a purely pseudoscalar mediator in the context of a 2HDM [44]. For our model, however, we have to generalise the results to arbitrary CP phases in the SM and dark sector. We use the opportunity to present intermediate calculational steps of the derivation not explicitly shown in ref. [44]. Our starting point is the calculation of the leading order effective vertices between the spin-0 mediator a and gluons, which we treat as background. To simplify the computation we employ the Fock-Schwinger gauge, in which the gluon field can directly be expressed in terms of the field strength [54,84]. The general amplitude can thus be written as where we have introduced the twist-2 gluon operator O G µν and a higher spin G a ρσ G aρσ (g αβ g µν − g αν g βµ ) .
(B.23) 12 The signs in front of O G µν and O G αβ differ from eq. (50) in ref. [54]. As we are not including the terms containing the twist-2 operator, this difference will not play a role.

JHEP06(2019)052
These operators do not contribute to G a ρσ G aρσ and give sub-leading SI interactions such that they are neglected in the present work. We therefore only have to consider the replacement G a αµ G a βν → 1 12 G a ρσ G aρσ (g αβ g µν − g αν g βµ ) , (B.24) for the spin-independent term. The spin-dependent term involving G a µν G aµν can be obtained straightforwardly by isolating the term in the trace containing the -tensor. Putting both contributions together, we then obtain (B.25) Since the momentum q corresponds to the second loop momentum of the full two-loop diagram it is not helpful to simply perform the loop integral over p, as one has to perform a second loop integral afterwards. It is more advantageous to use a Feynman parameter x instead: Shifting p → p + qx and performing the loop integral over p, we obtain finally One can proceed in a similar fashion for the remaining two terms in eq. (B.17), which contribute equally to the amplitude, such that the full amplitude reads where we have made the integral over the Feynman parameter symmetric under x → 1 − x. This result agrees with eq. (B.73) from ref. [44] for φ SM = π/2. Moreover, when performing a heavy quark expansion on eq. (B.29), we recover the effective Lagrangian from eq. (2.14).
It is now straightforward to perform the integration from the remaining triangle diagram of the full two-loop approach, visualised after the first arrow in figure 3, for which the amplitude can be written as Here we have introduced the shorthand notation

C Nuclear form factors
In this appendix we define the nuclear form factors required to calculate the effective interactions between DM and nucleons. For the spin-independent interactions we need the following nuclear form factors [56,85]: The second moments are calculated at the scale µ = m Z by using CTEQ PDFs [44,87].  13) where N refers to a change of nucleon momentum. This explicit dependence on the momentum transfer q µ arises from non-negligible π and η pole contributions. The corresponding expressions are given in eqs. (A30) and (A42) of ref. [88] and are implemented in DirectDM [60]. 13 We refer to ref. [86] for a discussion of the uncertainties of these form factors, in particular regarding the strange quark matrix element.

D UV-divergent loops in CP-violating Higgs portal
In this appendix we provide further details on the additional loop diagrams occurring in the CP-violating Higgs portal model, which are shown in figure 12. Like the diagrams considered in the main text, these diagrams contribute at the order 1/Λ 2 and induce spinindependent interactions for φ = π/2. The contribution of these diagrams to the triangle coefficients defined in eqs. (2.7) and (2.8) are given by where we have used the replacements from eq. (3.11). The loop integrals B 0 and B 1 are UV divergent and we replace the divergences by a logarithmic dependence on the new physics scale Λ from eq. (3.7) according to 1/ + ln(µ 2 /m 2 χ ) → ln(Λ 2 /m 2 χ ) [13]. We study the impact of this additional contribution in figure 13 for Λ = 1 TeV. We observe that, while the additional diagrams make the loop-contributions more important, the general conclusions drawn from figure 9 are not changed.