††institutetext:  Institute for Theoretical Particle Physics and Cosmology
(TTK), RWTH Aachen University,
D-52056 Aachen, Germany

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

Fatih Ertas  and Felix Kahlhoefer<EMAIL_ADDRESS><EMAIL_ADDRESS>
###### Abstract

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.

###### Keywords:

Mostly Weak Interactions: Beyond Standard Model; Astroparticles: Cosmology of
Theories beyond the SM

††preprint: TTK-19-08

##  1 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 Cushman:2013zza . Indeed, these experiments have become so
sensitive that they can constrain even DM models where standard spin-
independent and spin-dependent interactions are absent at leading order
Xia:2018qgs . As a result there has been a rapidly growing interest in DM
models with momentum- or velocity-dependent interactions, which can be
described by a general effective field theory (EFT) in the non-relativistic
limit Fitzpatrick:2012ix ; Anand:2013yka ; Catena:2014uqa ; Gresham:2014vja ;
Catena:2014epa ; Gluscevic:2015sqa ; Dent:2015zpa ; Kahlhoefer:2016eds ;
Bishara:2016hek ; Edwards:2018lsl . In these models it becomes essential to
include loop effects, which may reintroduce spin-independent interactions and
thereby substantially boost the expected event rates Haisch:2013uaa ;
Crivellin:2014qxa ; Crivellin:2014gpa ; DEramo:2016gos ; Bishara:2018vix .

Particular attention has been paid to models in which DM scattering is
mediated by a pseudoscalar exchange particle Freytsis:2010ne ; Dienes:2013xya
; Boehm:2014hva , motivated partially by the interesting implications for
collider Ipek:2014gua ; No:2015xqa ; Goncalves:2016iyg ; Bauer:2017ota ;
Bauer:2017fsw ; Pani:2017qyd ; Tunney:2017yfp ; Banerjee:2017wxi ;
Haisch:2018kqx ; Abe:2018bpo  and flavour Batell:2009di ; Freytsis:2009ct ;
Batell:2009jf ; Dolan:2014ska ; Berlin:2015wwa ; Dobrich:2018jyi  physics. At
leading order the resulting interactions are so strongly suppressed in the
non-relativistic limit that they are well below the “neutrino floor” which
indicates the ultimate reach of direct detection experiments Billard:2013qya .
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 Arcadi:2017wqi ; Sanderson:2018lmj ;
Li:2018qip ; Abe:2018emu . In fact, ref. Abe:2018emu  pointed out that for
this particular model even two-loop 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 Akerib:2018lyp  or XENONnT Aprile:2015uzo .
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 Kahlhoefer:2017umn ) loop effects
can be comparable to the leading-order 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 Beniwal:2015sdl ; Athron:2018hpc ; Abe:2019wku . 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.

##  2 Loop effects in direct detection

We investigate a simplified model of a Dirac fermion DM particle $\chi$
interacting with SM fermions $f$ through a general spin-0 mediator $a$ with
mass $m_{a}$ greater than the bottom-quark mass $m_{b}$:

$\displaystyle\mathcal{L}=g_{\chi}\,a\,\bar{\chi}\,(\cos\phi_{\chi}+i\gamma_{5}\sin\phi_{\chi})\,\chi+g_{\text{SM}}\sum_{f}\frac{m_{f}}{v}a\,\bar{f}\,(\cos\phi_{\text{SM}}+i\gamma_{5}\sin\phi_{\text{SM}})\,f\;,$
(1)

where $\phi_{\chi}$ and $\phi_{\text{SM}}$ are CP phases, $v\approx
246\,\mathrm{GeV}$ is the electroweak vacuum expectation value, $m_{f}$ are
the SM fermion masses and $g_{\chi}$ as well as $g_{\text{SM}}$ denote the
couplings of $a$ to DM and SM fermions, respectively. We have further assumed
Yukawa-like couplings in agreement with the hypothesis of minimal flavour
violation (MFV) DAmbrosio:2002vsn  such that flavour physics constraints on
the universal coupling $g_{\text{SM}}$ are weakened (see section 3.2).111In a
generic MFV scenario a slightly more general Lagrangian than eq. (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. For $\phi_{\chi}=\phi_{\text{SM}}=0$ we recover the well-
known simplified model of a scalar mediator, whereas for
$\phi_{\chi}=\phi_{\text{SM}}=\pi/2$ we obtain a CP-conserving theory with a
pseudoscalar mediator Abdallah:2015ter . 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 Arcadi:2017wqi ; Sanderson:2018lmj ;
Abe:2018emu . 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. (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 $\lambda_{ah}$ and supplement eq. (1) by the interaction term

$\displaystyle\mathcal{L}^{\text{Higgs}}_{\text{int}}=\frac{1}{2}\lambda_{ah}vha^{2}\;.$
(2)

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.

###  2.1 Low-energy effective Lagrangian

Figure 1: Tree-level, Higgs-induced triangle as well as box diagram
contribution to the cross section relevant for direct searches of DM. All
Feynman diagrams are drawn with TikZ-Feynman Ellis:2016jkw .

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.222Note
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. Starting
with the tree-level exchange of $a$ illustrated in the left panel of figure 1,
we obtain

$\displaystyle\mathcal{L}^{\text{SI}}_{\text{tree}}$
$\displaystyle=\sum_{q=\text{all}}m_{q}\,\mathcal{C}^{\text{tree}}\big{(}\cos(\phi_{\chi})\cos(\phi_{\text{SM}})\,\bar{\chi}\chi+\sin(\phi_{\chi})\cos(\phi_{\text{SM}})\,\bar{\chi}i\gamma_{5}\chi\big{)}\,\bar{q}q\;,$
(3) $\displaystyle\mathcal{L}^{\text{SD}}_{\text{tree}}$
$\displaystyle=\sum_{q=\text{all}}m_{q}\,\mathcal{C}^{\text{tree}}\big{(}\cos(\phi_{\chi})\sin(\phi_{\text{SM}})\,\bar{\chi}\chi+\sin(\phi_{\chi})\sin(\phi_{\text{SM}})\,\bar{\chi}i\gamma_{5}\chi\big{)}\,\bar{q}i\gamma_{5}q\;,$
(4)

where the sum runs over all quark species. Here we have defined the tree-level
coefficient

$\displaystyle\mathcal{C}^{\text{tree}}=\frac{g_{\chi}\,g_{\text{SM}}}{v\,m_{a}^{2}}\;,$
(5)

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

$\displaystyle\mathcal{L}^{\text{SI}}_{\text{triangle}}=\sum_{q=\text{all}}\frac{m_{q}\lambda_{ah}}{m_{h}^{2}}\left(\mathcal{C}_{S}^{\text{triangle}}\,\bar{\chi}\chi\,\bar{q}q+\mathcal{C}_{PS}^{\text{triangle}}\,\bar{\chi}i\gamma_{5}\chi\,\bar{q}q\right)\;,$
(6)

where the sum again runs over all quarks. We have further introduced the
triangle coefficients

$\displaystyle\mathcal{C}_{S}^{\text{triangle}}$
$\displaystyle=\frac{g_{\chi}^{2}}{(4\pi)^{2}}\,m_{\chi}\left[(1+\cos(2\phi_{\chi}))\,C_{0}(m_{\chi}^{2},\,m_{a}^{2},\,m_{\chi}^{2})+C_{2}(m_{\chi}^{2},\,m_{a}^{2},\,m_{\chi}^{2})\right]\;,$
(7) $\displaystyle\mathcal{C}_{PS}^{\text{triangle}}$
$\displaystyle=\frac{g_{\chi}^{2}}{(4\pi)^{2}}\,m_{\chi}\sin(2\phi_{\chi})\,C_{0}(m_{\chi}^{2},\,m_{a}^{2},\,m_{\chi}^{2})\;,$
(8)

in terms of the loop functions
$C_{0}(m_{\chi}^{2},\,m_{a}^{2},\,m_{\chi}^{2})$ and
$C_{2}(m_{\chi}^{2},\,m_{a}^{2},\,m_{\chi}^{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 Abe:2018emu , and obtain

$\displaystyle\begin{split}\mathcal{L}^{\text{SI}}_{\text{box}}=&\sum_{q=u,d,s}\left(m_{q}\,\mathcal{C}^{\text{box}}_{1,q}\,\bar{\chi}\chi\,\bar{q}q+m_{q}\,\mathcal{C}^{\text{box}}_{2,q}\,\bar{\chi}i\gamma_{5}\chi\,\bar{q}q\right)\\\
&+\sum_{q=u,d,s,c,b}\Big{(}\mathcal{C}^{\text{box}}_{5,q}\,\bar{\chi}i\partial^{\mu}\gamma^{\nu}\chi\,\mathcal{O}^{q}_{\mu\nu}+\mathcal{C}^{\text{box}}_{6,q}\,\bar{\chi}i\partial^{\mu}i\partial^{\nu}\chi\,\mathcal{O}^{q}_{\mu\nu}+\mathcal{C}^{\text{box}}_{7,q}\,\bar{\chi}i\gamma_{5}i\partial^{\mu}i\partial^{\nu}\chi\,\mathcal{O}^{q}_{\mu\nu}\Big{)}\;,\end{split}$
(9) $\displaystyle\mathcal{L}^{\text{SD}}_{\text{box}}=$
$\displaystyle\sum_{q=u,d,s}\left(m_{q}\,\mathcal{C}^{\text{box}}_{3,q}\,\bar{\chi}\chi\,\bar{q}i\gamma_{5}q+m_{q}\,\mathcal{C}^{\text{box}}_{4,q}\,\bar{\chi}i\gamma_{5}\chi\,\bar{q}i\gamma_{5}q\right)\;.$
(10)

Computational details and the expressions of the different box diagram
coefficients $\mathcal{C}^{\text{box}}_{i,q}$ are given in appendix A.2. Note
that all of these coefficients share a common factor of
$g^{2}_{\chi}\,g^{2}_{\text{SM}}\,m_{q}^{2}/v^{2}$, which also constitutes the
only quark dependence. In eq. (9) we have also introduced the twist-2 quark
operator

$\displaystyle\mathcal{O}^{q}_{\mu\nu}=\bar{q}\,\left(\frac{i\partial^{\mu}\gamma^{\nu}+i\partial^{\nu}\gamma^{\mu}}{2}-\frac{1}{4}g^{\mu\nu}i\not{\partial}\right)\,q\;.$
(11)

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. (9) Hisano:2010ct ; Hisano:2015bma . However, none of the heavy
quarks have been included in the remaining terms of eqs. (9) and (10), because
they require a different treatment, which will be discussed next.

###  2.2 Effective description of two-loop processes

Figure 2: Two-loop processes for the evaluation of the heavy quark ($Q=c,b,t$)
contributions to effective DM-gluon interactions.

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 Shifman:1978zn

$\displaystyle m_{Q}\bar{Q}Q$
$\displaystyle\rightarrow-\frac{\alpha_{s}}{12\pi}G^{a}_{\mu\nu}G^{a\mu\nu}\;,$
(12) $\displaystyle m_{Q}\bar{Q}i\gamma_{5}Q$
$\displaystyle\rightarrow\frac{\alpha_{s}}{8\pi}G^{a}_{\mu\nu}\widetilde{G}^{a\mu\nu}\;,$
(13)

where $G^{a\mu\nu}$ is the gluon field strength tensor and
${\widetilde{G}}^{a\mu\nu}=\frac{1}{2}\epsilon^{\mu\nu\alpha\beta}G^{a}_{\alpha\beta}$
with the convention $\epsilon^{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. (12) would
neglect the contribution from the diagram on the right. For $m_{Q}\ll
m_{a},\,m_{\chi}$ the two-loop computation hence cannot be simplified in this
way without introducing potentially large errors Abe:2018emu . In the opposite
case of $m_{Q}\gg m_{a},\,m_{\chi}$ it was argued in ref. Abe:2018emu  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.

Figure 3: Illustration of the decomposition of the two-loop process for
$m_{Q}\gg m_{a},\,m_{\chi}$. 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.

\+ $\rightarrow$ $\rightarrow$

As we are now going to demonstrate, however, for $m_{Q}\gg m_{a},\,m_{\chi}$
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}\ll 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 Patel:2015tea  for the
evaluation and expansion of the loop computations. This then yields the
following leading order effective Lagrangian coupling $a$ to gluons

$\displaystyle\mathcal{L}_{\text{eff}}^{\text{aaG}}$
$\displaystyle=\frac{1}{2}\,d^{\text{\,eff}}_{G}\,aa\,\frac{\alpha_{s}}{12\pi}\,G^{a}_{\mu\nu}G^{a\mu\nu}+\frac{1}{2}\,d^{\text{\,eff}}_{\widetilde{G}}\,aa\,\frac{\alpha_{s}}{8\pi}G^{a}_{\mu\nu}\widetilde{G}^{a\mu\nu}\;.$
(14)

Here we have included a symmetry factor of 1/2 and defined333Note that
$d^{\text{\,eff}}_{G}$ vanishes for certain values of $\phi_{\text{SM}}$ such
that one would need to include higher orders. However, these specific cases
are not of interest in the present work. While
$d^{\text{\,eff}}_{\widetilde{G}}$ also vanishes for specific values of
$\phi_{\text{SM}}$, the same is true for the full expression
$d^{\text{\,full}}_{\widetilde{G}}$, see eq. (90) in appendix B.2, i.e. this
is not a result of the heavy quark expansion.

$\displaystyle d^{\text{\,eff}}_{G}=\frac{g_{\text{SM}}^{2}}{v^{2}}$
$\displaystyle\big{(}\sin^{2}(\phi_{\text{SM}})-\cos^{2}(\phi_{\text{SM}})\big{)}\;,\qquad
d^{\text{\,eff}}_{\widetilde{G}}=\frac{g_{\text{SM}}^{2}}{v^{2}}\sin(2\phi_{\text{SM}})\;,$
(15)

which are both independent of the top-quark mass. Now performing the second
step visualised in figure 3, we obtain for the effective two-loop approach

$\displaystyle\mathcal{L}^{\text{SI}}_{\text{2-Loop}}$
$\displaystyle=\left(\mathcal{C}^{\text{eff}}_{G,S}\,\bar{\chi}\chi+\mathcal{C}^{\text{eff}}_{G,PS}\,\bar{\chi}i\gamma_{5}\chi\,\right)\frac{-\alpha_{s}}{12\pi}\,G^{a}_{\mu\nu}G^{a\mu\nu}\;,$
(16) $\displaystyle\mathcal{L}^{\text{SD}}_{\text{2-Loop}}$
$\displaystyle=\left(\mathcal{C}^{\text{eff}}_{\widetilde{G},S}\,\bar{\chi}\chi+\mathcal{C}^{\text{eff}}_{\widetilde{G},PS}\,\bar{\chi}i\gamma_{5}\chi\,\right)\frac{\alpha_{s}}{8\pi}\,G^{a}_{\mu\nu}{\widetilde{G}}^{a\mu\nu}\;,$
(17)

where the effective two-loop coefficients read

$\displaystyle\mathcal{C}^{\text{eff}}_{G,S}$
$\displaystyle=d^{\text{\,eff}}_{G}\,\mathcal{C}_{S}^{\text{triangle}}\;,\qquad$
$\displaystyle\mathcal{C}^{\text{eff}}_{G,PS}$
$\displaystyle=d^{\text{\,eff}}_{G}\,\mathcal{C}_{PS}^{\text{triangle}}\;,$
(18) $\displaystyle\mathcal{C}^{\text{eff}}_{G,S}$
$\displaystyle=-d^{\text{\,eff}}_{\widetilde{G}}\,\mathcal{C}_{S}^{\text{triangle}}\;,\qquad$
$\displaystyle\mathcal{C}^{\text{eff}}_{\widetilde{G},PS}$
$\displaystyle=-d^{\text{\,eff}}_{\widetilde{G}}\,\mathcal{C}_{PS}^{\text{triangle}}\;.$
(19)

An analogous calculation for the bottom and charm quark only gives a useful
approximation if $m_{a}\ll 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.444This conclusion
could change for example in models with extended Higgs sectors, where
couplings to down-type quarks may receive a substantial enhancement. 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.

Figure 4: Comparison of $|\mathcal{C}_{G,S}|$ in the effective approach for
the top quark (red), the two-loop contribution of the top quark (dotted green)
and the two-loop result including all heavy quarks (dashed grey) as a function
of $m_{a}$ (left panel) and $m_{\chi}$ (right panel). For both plots we set
$\phi_{\chi}=0$ and $\phi_{\text{SM}}=\pi/2$.

This is illustrated in figure 4, where we plot the absolute value of the
coefficient $\mathcal{C}_{G,S}$ as a function of the mediator mass (left
panel) and of the DM mass $m_{\chi}$ (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}\ll m_{t}$
across the whole range of DM masses. Including also bottom and charm quark in
the two-loop calculation has only slight influences for small values of
$m_{a}$ (dashed grey). Similar results can be obtained for the other
coefficients.555For specific parameter points cancellations might occur within
the coefficients $\mathcal{C}_{G}$ and $\mathcal{C}_{\widetilde{G}}$ like in
$d_{G}^{\text{\,eff}}$ for $\phi_{\text{SM}}\approx\pi/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. We
conclude that it is possible to simplify the full two-loop calculation in the
case of $m_{Q}\gg m_{a},\,m_{\chi}$, 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.

###  2.3 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:

$\displaystyle\mathcal{L}_{\chi N}^{\text{eff}}$
$\displaystyle=\left(\mathcal{C}^{\text{SI}}_{\text{eff},N}\,\bar{\chi}\chi+\mathcal{C}^{\text{SI,CPV}}_{\text{eff},N}\,\bar{\chi}i\gamma_{5}\chi\right)\bar{N}N+\left(\mathcal{C}^{\text{SD,CPV}}_{\text{eff},N}\,\bar{\chi}\chi+\mathcal{C}^{\text{SD}}_{\text{eff},N}\,\bar{\chi}i\gamma_{5}\chi\right)\bar{N}i\gamma_{5}N\;,$
(20)

where $N=p,n$ is a nucleon field and ‘CPV’ indicates terms that only arise
when CP is violated. The coefficients $C_{\text{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_{\text{eff},p}\neq
C_{\text{eff},n}$.

For the SI coefficients, we find

$\displaystyle\mathcal{C}^{\text{SI}}_{\text{eff},N}$
$\displaystyle=\sum_{q=u,d,s}m_{N}f^{N}_{q}\left(\cos(\phi_{\chi})\cos(\phi_{\text{SM}})\,\mathcal{C}^{\text{tree}}+\frac{\lambda_{ah}}{m_{h}^{2}}\,\mathcal{C}^{\text{triangle}}_{S}+\mathcal{C}^{\text{box}}_{1,q}\right)$
$\displaystyle\hskip
14.22636pt+3\cdot\frac{2}{27}m_{N}f^{N}_{G}\left(\cos(\phi_{\chi})\cos(\phi_{\text{SM}})\,\mathcal{C}^{\text{tree}}+\frac{\lambda_{ah}}{m_{h}^{2}}\,\mathcal{C}^{\text{triangle}}_{S}\right)$
(21) $\displaystyle\hskip
14.22636pt+\sum_{q=u,d,s,c,b}\frac{3}{4}m_{N}m_{\chi}\,\Big{(}q^{N}(2)+\bar{q}^{N}(2)\Big{)}\Big{(}\mathcal{C}^{\text{box}}_{5,q}+m_{\chi}\,\mathcal{C}^{\text{box}}_{6,q}\Big{)}+\frac{2}{27}m_{N}f^{N}_{G}\,\mathcal{C}^{\text{eff}}_{G,S}\;,$

as well as

$\displaystyle\mathcal{C}^{\text{SI,CPV}}_{\text{eff},N}$
$\displaystyle=\sum_{q=u,d,s}m_{N}f^{N}_{q}\left(\sin(\phi_{\chi})\cos(\phi_{\text{SM}})\,\mathcal{C}^{\text{tree}}+\frac{\lambda_{ah}}{m_{h}^{2}}\,\mathcal{C}^{\text{triangle}}_{PS}+\mathcal{C}^{\text{box}}_{2,q}\right)$
$\displaystyle\hskip
14.22636pt+3\cdot\frac{2}{27}m_{N}f^{N}_{G}\left(\sin(\phi_{\chi})\cos(\phi_{\text{SM}})\,\mathcal{C}^{\text{tree}}+\frac{\lambda_{ah}}{m_{h}^{2}}\,\mathcal{C}^{\text{triangle}}_{PS}\right)$
(22) $\displaystyle\hskip
14.22636pt+\sum_{q=u,d,s,c,b}\frac{3}{4}m_{N}m^{2}_{\chi}\,\Big{(}q^{N}(2)+\bar{q}^{N}(2)\Big{)}\,\mathcal{C}^{\text{box}}_{7,q}+\frac{2}{27}m_{N}f^{N}_{G}\,\mathcal{C}^{\text{eff}}_{G,PS}\;,$

where the nuclear form factors $f^{N}_{q,G}$, $q^{N}(2)$ and $\bar{q}^{N}(2)$
are defined in appendix C.

Likewise, we obtain for the SD coefficients

$\displaystyle\begin{split}\mathcal{C}^{\text{SD,CPV}}_{\text{eff},N}&=\sum_{q=u,d,s}F_{P}^{q/N}(q^{2})\,\Big{(}\cos(\phi_{\chi})\sin(\phi_{\text{SM}})\,\mathcal{C}^{\text{tree}}+\mathcal{C}^{\text{box}}_{3,q}\Big{)}\\\
&\hskip
14.22636pt+\,F^{N}_{\widetilde{G}}(q^{2})\,\Big{(}3\cos(\phi_{\chi})\sin(\phi_{\text{SM}})\,\mathcal{C}^{\text{tree}}+\,\mathcal{C}^{\text{eff}}_{\widetilde{G},S}\Big{)}\;,\end{split}$
(23)

and

$\displaystyle\begin{split}\mathcal{C}^{\text{SD}}_{\text{eff},N}&=\sum_{q=u,d,s}F_{P}^{q/N}(q^{2})\,\Big{(}\sin(\phi_{\chi})\sin(\phi_{\text{SM}})\,\mathcal{C}^{\text{tree}}+\mathcal{C}^{\text{box}}_{4,q}\Big{)}\\\
&\hskip
14.22636pt+\,F^{N}_{\widetilde{G}}(q^{2})\,\Big{(}3\sin(\phi_{\chi})\sin(\phi_{\text{SM}})\,\mathcal{C}^{\text{tree}}+\,\mathcal{C}^{\text{eff}}_{\widetilde{G},PS}\Big{)}\;.\end{split}$
(24)

The form factors $F_{P}^{q/N}(q^{2})$ and $F^{N}_{\widetilde{G}}(q^{2})$ are
given in appendix C. Because of non-negligible contributions from the $\pi$
and $\eta$ pole, these form factors depend on the momentum exchange $q^{\mu}$
between DM and nucleons.

In the non-relativistic limit the effective Lagrangian from eq. (20) can be
matched onto a basis of effective operators:

$\mathcal{L}_{\chi N}^{\text{eff}}\to\sum_{i}c_{i}^{N}\mathcal{O}^{N}_{i}\;,$
(25)

where the operators $\mathcal{O}^{N}_{i}$ depend only on the spins
$\vec{S}_{\chi}$ and $\vec{S}_{N}$ of DM and the nucleon, respectively, as
well as on the momentum transfer $\vec{q}$ and the DM-nucleon relative
velocity $\vec{v}$ Fitzpatrick:2012ix ; Anand:2013yka ; Fan:2010gt . For the
model that we consider, only four different operators are generated, namely

$\displaystyle\begin{split}\mathcal{O}^{N}_{1}&=1\,,\\\
\mathcal{O}^{N}_{6}&=(\vec{S}_{\chi}\cdot\frac{\vec{q}}{m_{N}})(\vec{S}_{N}\cdot\frac{\vec{q}}{m_{N}})\,,\\\
\mathcal{O}^{N}_{10}&=i(\vec{S}_{N}\cdot\frac{\vec{q}}{m_{N}})\,,\\\
\mathcal{O}^{N}_{11}&=i(\vec{S}_{\chi}\cdot\frac{\vec{q}}{m_{N}})\;.\end{split}$
(26)

The corresponding coefficients can be directly read off from
$\mathcal{L}_{\chi N}^{\text{eff}}$ Anand:2013yka :

$c_{1}^{N}=C_{\text{eff},N}^{\text{SI}}\,,\qquad
c_{6}^{N}=\frac{m_{N}}{m_{\chi}}C_{\text{eff},N}^{\text{SD}}\,,\qquad
c_{10}^{N}=C_{\text{eff},N}^{\text{SD,CPV}}\,,\qquad
c_{11}^{N}=-\frac{m_{N}}{m_{\chi}}C_{\text{eff},N}^{\text{SI,CPV}}\;.$ (27)

Note that like the form factors $F_{P}^{q/N}(q^{2})$ and
$F^{N}_{\widetilde{G}}(q^{2})$ the coefficients $c_{6}^{N}$ and $c_{10}^{N}$
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. (1).

##  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

$\sigma^{\text{SI}}_{N}=\frac{\mu_{N}^{2}\,|c_{1}^{N}|^{2}}{\pi}\;,$ (28)

where $\mu_{N}=m_{\chi}m_{N}/(m_{\chi}+m_{N})$ is the DM-nucleon reduced mass.
For $c_{1}^{p}\approx c_{1}^{n}$ the differential event rate with respect to
recoil energy $E_{\mathrm{R}}$ is then simply given by

$\displaystyle\frac{\mathrm{d}R}{\mathrm{d}E_{\mathrm{R}}}=\frac{\rho_{0}\,\sigma^{\text{SI}}_{p}\,A^{2}\,F^{2}(E_{\text{R}})}{2\,\mu_{p}^{2}\,m_{\chi}}g(v_{\text{min}})\;,$
(29)

where $\rho_{0}$ is the local DM density, $A$ is the mass number of the target
nucleus and $F^{2}(E_{\text{R}})$ denotes the nuclear form factor. The factor
$g(v_{\text{min}})=\int_{v_{\text{min}}}f(v)/v\,\mathrm{d}v$ denotes the
velocity integral as a function of the minimum velocity
$v_{\text{min}}(E_{\text{R}})=\sqrt{m_{A}E_{\text{R}}/(2\,\mu^{2})}$ with
$m_{A}$ being the mass of the target nucleus and $\mu$ 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 $\sigma^{\text{SI}}_{p}$ as a
function of $m_{\chi}$.

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. Fitzpatrick:2012ix ;
Anand:2013yka ; Kahlhoefer:2016eds . 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 Workgroup:2017lvb ;
Athron:2018hpc , which includes an extensive database of existing and planned
direct detection experiments. Furthermore, DDCalc contains an automated
interface with DirectDM Bishara:2017nnn , which we use for the matching of the
spin-dependent coefficients in eq. (27) and the evaluation of the
corresponding nuclear form factors. We can therefore simply pass the
coefficients $C_{\text{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 Aprile:2018dbl  and the regions that predict at least 5 events in the
next-generation LZ experiment Akerib:2018lyp .666This number of events
corresponds approximately to the median expected sensitivity using a cut-and-
count analysis with a background expectation of 6.49 events. A better
sensitivity may be achieved by exploiting differences in the differential
distributions between signal and background. Similar exclusion limits are
obtained from the Panda-X Cui:2017nnn ; Xia:2018qgs  and LUX Akerib:2016vxi ;
Akerib:2017kat  experiments, while comparable sensitivities are expected for
the XENONnT experiment Aprile:2015uzo .

###  3.1 General CP phases

Figure 5:  Direct detection constraints as a function of the CP-violating
phases $\phi_{\chi}$ and $\phi_{\text{SM}}$. For both panels we have fixed
$g_{\chi}=g_{\text{SM}}=1$ and $\lambda_{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.

We first consider in figure 5 the most general case, in which $\phi_{\chi}$
and $\phi_{\text{SM}}$ can take arbitrary values between 0 (corresponding to
purely scalar couplings) and $\pi/2$ (corresponding to purely pseudoscalar
couplings). For the purpose of this figure we have fixed
$g_{\text{SM}}=g_{\chi}=1$ and $\lambda_{ah}=0$ and consider two different
combinations of $m_{\chi}$ and $m_{a}$ in the two panels. 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 $\phi_{\chi}$, $\phi_{\text{SM}}\ll\pi/2$ the tree-level exchange of $a$
dominates the spin-independent coefficient
$\mathcal{C}^{\text{SI}}_{\text{eff},N}$ from eq. (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_{\text{SM}}$ to
be very small Kaplinghat:2013yxa ; Kahlhoefer:2017umn . As the two phases
approach $\pi/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 $\phi_{\chi}=0$ and $\phi_{\text{SM}}=\pi/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 $\mathcal{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 $\mathcal{O}(1)$ couplings Dienes:2013xya . 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 $\phi_{\chi}\approx 0$ and
$10^{-4}\lesssim\pi/2-\phi_{\text{SM}}\lesssim 10^{-3}$, where the total event
rate is _smaller_ than the one predicted at tree-level due to the destructive
interference between spin-independent 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 $\phi_{\chi}=\pi/2$ and $\phi_{\text{SM}}=0$,
i.e. the bottom-right corner of figure 5, the tree-level contribution to spin-
independent scattering maps onto the non-relativistic operator
$\mathcal{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_{N}^{2}/m_{\chi}^{2}$. We will therefore study the influence of purely
spin-independent contributions emerging at loop-level in the context of the
CP-violating 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 Abe:2018emu , 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.

###  3.2 Maximal CP violation

Let us take a closer look at the case $\phi_{\chi}=0$ and
$\phi_{\text{SM}}=\pi/2$, corresponding to the top-left 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 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_{\text{SM}}$, $g_{\chi}$ and $\lambda_{ah}$.

For given values of $m_{a}$, $m_{\chi}$ and $g_{\text{SM}}$ we can fix
$g_{\chi}$ by the requirement that the observed DM relic abundance can be
explained in terms of thermal freeze-out via the annihilation processes
$\chi\bar{\chi}\to q\bar{q}$ and $\chi\bar{\chi}\to aa$. If the latter process
is kinematically allowed, i.e. for $m_{a}<m_{\chi}$, it will typically give
the dominant contribution for $g_{\text{SM}}\ll 1$, such that the required
value for $g_{\chi}$ becomes independent of $g_{\text{SM}}$. In this limit, we
find $g_{\chi}\propto m_{\chi}^{1/2}$ with $g_{\chi}=1$ for $m_{\chi}\approx
500\,\mathrm{GeV}$. For larger $g_{\text{SM}}$ the calculation becomes more
involved and we use micrOmegas_v5.0.6 Belanger:2018mqt  to determine the
required value for $g_{\chi}$ 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}\lesssim
m_{B}\approx 5.2\,\mathrm{GeV}$, these constraints are very strong and
effectively exclude the possibility of obtaining observable direct detection
signatures Dolan:2014ska . However, almost all of these constraints disappear
for larger values of $m_{a}$. Bounds from radiative $\Upsilon$ decays
Lees:2012iw ; Lees:2011wb  extend to slightly larger masses, but also
disappear for $m_{a}\gtrsim 7\,\mathrm{GeV}$. Provided the pseudoscalar
couples also to leptons (with coupling strength $g_{\text{SM}}\,m_{\ell}/v$),
another important constraint arises from $B_{s}\to\mu^{+}\mu^{-}$, which can
arise from loop-induced flavour-changing interactions with an off-shell
mediator. The resulting branching ratio is given by Altmannshofer:2011gn ;
Dolan:2014ska

$\displaystyle\frac{\text{BR}(B_{s}\rightarrow\mu^{+}\mu^{-})_{\text{NP}}}{\text{BR}(B_{s}\rightarrow\mu^{+}\mu^{-})_{\text{SM}}}$
$\displaystyle\simeq\frac{g_{\text{SM}}^{4}\,m_{t}^{4}\,m_{B_{s}}^{4}}{16m_{W}^{4}\,\sin(\theta_{W})^{4}|C_{10}^{\text{SM}}|^{2}\,\left((m_{B_{s}}^{2}-m_{a}^{2})^{2}+\Gamma_{a}^{2}\,m_{a}^{2}\right)}\log^{2}\left(\frac{\Lambda^{2}}{m_{t}^{2}}\right)\;,$
(30)

where $C_{10}^{\text{SM}}=-4.103$ and $\Lambda$ is the scale of new physics
(such as additional charged Higgs bosons needed in a gauge-invariant UV
completion). For $\Lambda=1\,\mathrm{TeV}$ and assuming $m_{a}\gg m_{B_{s}}$,
this expression simplifies to

$\displaystyle\frac{\text{BR}(B_{s}\rightarrow\mu^{+}\mu^{-})_{\text{NP}}}{\text{BR}(B_{s}\rightarrow\mu^{+}\mu^{-})_{\text{SM}}}\approx\left(\frac{7.9\,\mathrm{GeV}\,g_{\text{SM}}}{m_{a}}\right)^{4}\;.$
(31)

The branching ratio of $B_{s}\to\mu^{+}\mu^{-}$ has been measured with a
precision of $20\%$ Aaij:2017vad  and is found to be in agreement with the SM
prediction Bobeth:2013uxa . To obtain an approximate bound on $g_{\text{SM}}$
we therefore require the new-physics contribution not to exceed 40% of the SM
value. This gives

$g_{\text{SM}}\lesssim 1.0\left(\frac{m_{a}}{10\,\mathrm{GeV}}\right)\;.$ (32)

In other words, even for spin-0 bosons as light as $10\,\mathrm{GeV}$ the
coupling strength $g_{\text{SM}}$ can be of order unity.777For
$\phi_{\text{SM}}$ different from $\pi/2$ there would be additional
constraints from observables sensitive to CP-violation, in particular electric
dipole moments of leptons Chen:2015vqy ; Marciano:2016yhf  and nuclei
Mantry:2014zsa . However, for $\phi_{\text{SM}}\approx\pi/2$ the spin-0
mediator behaves like a pure pseudoscalar in all observables involving only SM
particles, such that these constraints are absent.

The situation is quite different for the coupling $\lambda_{ah}$ between $a$
and the SM Higgs boson. This coupling induces the decay $h\to aa$ with partial
width Beniwal:2015sdl

$\Gamma_{h\to
aa}=\frac{\lambda_{ah}^{2}\,v^{2}}{32\pi\,m_{h}}\left(1-\frac{4\,m_{a}^{2}}{m_{h}^{2}}\right)^{1/2}\;.$
(33)

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 Haisch:2018kqx ,
at present the strongest constraints come from a global fit of the measured
properties of the SM-like Higgs boson at ATLAS and CMS Khachatryan:2016vau .
These fits imply $\text{BR}(h\to aa)<0.34$, corresponding to $\Gamma_{h\to
aa}\lesssim 2\,\mathrm{MeV}$, when simultaneously allowing for modifications
of the Higgs boson production cross section, or $\text{BR}(h\to aa)<0.13$,
corresponding to $\Gamma_{h\to aa}\lesssim 0.6\,\mathrm{MeV}$, when assuming
the production cross section to be given by the SM prediction. For $m_{a}\ll
m_{h}/2$, these bounds translate to $\lambda_{ah}\lesssim 0.02$ and
$\lambda_{ah}\lesssim 0.01$, respectively. We will conservatively show the
weaker bound in the following.

Figure 6:  Constraints on $g_{\text{SM}}$ (left) and $\lambda_{ah}$ (right) as
a function of $m_{\chi}$ in a model with maximal CP violation
($\phi_{\chi}=0$, $\phi_{\text{SM}}=\pi/2$). At each point the coupling
$g_{\chi}$ is fixed in such a way that the observed DM relic abundance is
reproduced. The dotted lines in the left panel indicate the ratio of the
predicted number of events in LZ from loop-induced spin-independent
interactions and from tree-level momentum-suppressed interactions. In the
right panel, tree-level interactions are absent.

Figure 6 summarises the constraints on $g_{\text{SM}}$ (left) and
$\lambda_{ah}$ (right) as a function of $m_{\chi}$. At each point in the two
plots $g_{\chi}$ is determined by the relic density requirement and we have
set $m_{a}=15\,\mathrm{GeV}$. Again the solid blue region is excluded by
XENON1T and the parameter points for which 5 events are predicted in LZ are
indicated by the dashed green line. Dotted black lines in the left panel
indicate the ratio of loop-induced spin-independent 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_{\text{SM}}$ and with increasing $m_{\chi}$, corresponding to
increasing $g_{\chi}$. The kinks for $m_{\chi}\approx 175\,\mathrm{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_{\chi}$. In the right panel, tree-level
interactions are absent (since we have set $g_{\text{SM}}=0$), so direct
detection constraints arise exclusively from loop-induced interactions.

Since direct constraints on $g_{\text{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 $\lambda_{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
$\phi_{\text{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_{\chi}=200\,\mathrm{GeV}$, $m_{a}=15\,\mathrm{GeV}$ and
$g_{\text{SM}}=0.6$, corresponding to $g_{\chi}=0.6$. The tree-level
interactions are momentum-suppressed and therefore vanish in the limit
$E_{\mathrm{R}}\to 0$, leading to a maximum around $E_{\mathrm{R}}\sim
30\,\mathrm{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
$\phi_{\chi}=\pi/2$, $\phi_{\text{SM}}=0$, in which case the tree-level
process is coherently enhanced but suppressed by a factor $m_{N}/m_{\chi}$ in
$c_{11}$, see eq. (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.

Figure 7:  Predicted differential event rate in LZ for a specific parameter
point in the model with maximal CP violation ($\phi_{\chi}=0$,
$\phi_{\text{SM}}=\pi/2$) consistent with all current constraints.

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 ($\phi_{\chi}=0$,
$\phi_{\text{SM}}=\pi/2$), the right panel corresponds to the pure
pseudoscalar case ($\phi_{\chi}=\phi_{\text{SM}}=\pi/2$). In both cases we fix
$g_{\chi}$ by the relic density requirement and set $g_{\text{SM}}=0.6$,
consistent with the bounds discussed above (which are independent of
$\phi_{\chi}$). 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. Arcadi:2017wqi ; Sanderson:2018lmj
.888Here we have used the coefficient $C_{S,q}$ from ref. Arcadi:2017wqi  for
the top quark and have fixed the overall sign following ref. Abe:2018emu  As
pointed out previously Abe:2018emu , this approach leads to a vast
overestimation of the loop contribution.

Figure 8: Comparison of the effective approach and the full two-loop result
for two benchmark points in $m_{\chi}-\sigma^{\text{SI}}_{p}\,-$ plane. In the
left panel we consider maximal CP violation with $\phi_{\text{SM}}=\pi/2$ and
$\phi_{\chi}=0$, whereas in the right panel we fix
$\phi_{\text{SM}}=\phi_{\chi}=\pi/2$, i.e. pure pseudoscalar phases, and also
show the curve corresponding to previous calculational approaches of the two-
loop diagram. In both panels $\lambda_{ah}$ is set to zero and $g_{\chi}$ 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.

###  3.3 CP-violating Higgs portal

As a final example for the importance of loop-effects we consider the
fermionic Higgs portal model Beniwal:2015sdl ; Athron:2018hpc ; Abe:2019wku :

$\mathcal{L}=\mathcal{L}_{\text{SM}}+\overline{\chi}(i\not{\partial}-\mu)\chi-\frac{\lambda_{h\chi}}{\Lambda}\left(\cos\psi\,\overline{\chi}\chi+\sin\psi\,\overline{\chi}i\gamma_{5}\chi\right)H^{\dagger}H\;,$
(34)

where $H$ denotes the SM Higgs doublet and $\Lambda$ parametrises the unknown
scale of new physics. At first sight, this Lagrangian bears little resemblance
to the simplified model discussed so far. After electroweak symmetry breaking,
however, the following interactions are generated:999We omit the interaction
of DM with two Higgs bosons, which plays no role for direct detection. This
interaction is however included in the relic density calculation presented
below.

$\mathcal{L}\supset-\lambda\,v\,h^{3}-\frac{h}{v}\sum_{q}m_{q}\overline{q}q-\frac{\lambda_{h\chi}\,v}{\Lambda}h\left(\cos\phi\,\overline{\chi}\chi+\sin\phi\,\overline{\chi}i\gamma_{5}\chi\right)\,,$
(35)

where $\lambda$ denotes the quartic Higgs self-coupling and

$\cos\phi=\frac{\mu}{m_{\chi}}\left(\cos\psi+\frac{1}{2}\frac{\lambda_{h\chi}}{\Lambda}\frac{v^{2}}{\mu_{\chi}}\right)\;,$
(36)

with

$\displaystyle m_{\chi}$
$\displaystyle=\sqrt{\left(\mu+\frac{1}{2}\frac{\lambda_{h\chi}}{\Lambda}v^{2}\cos\psi\right)^{2}+\left(\frac{1}{2}\frac{\lambda_{h\chi}}{\Lambda}v^{2}\sin\psi\right)^{2}}\,.$
(37)

We can therefore directly apply all the results from section 2 with the
replacements

$\displaystyle m_{a}=m_{h},\quad
g_{\chi}=\frac{\lambda_{h\chi}v}{\Lambda},\quad
g_{\text{SM}}=1,\quad\phi_{\text{SM}}=0,\quad\phi_{\chi}=\phi,\quad\lambda_{ah}=-6\lambda=-3\frac{m_{h}^{2}}{v^{2}}\;.$
(38)

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_{\chi}$, $\lambda_{h\chi}/\Lambda$ and
$\phi$.

For $\phi\neq 0$ the model violates CP and spin-independent scattering is
suppressed proportional to $\cos^{2}\phi$. As $\phi$ approaches $\pi/2$, loop
effects are therefore expected to become increasingly important. We confirm
this expectation in figure 9, which shows constraints on
$\lambda_{h\chi}/\Lambda$ as a function of $m_{\chi}$. 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.101010We note that our effective description of top-
quark loops overestimates the contribution to the Wilson coefficient for spin-
independent scattering by up to a factor of 3 compared to the full two-loop
result. However, by far the dominant contribution to this coefficient arises
from triangle diagrams, making the difference between the effective
description and the full two-loop calculation irrelevant.

Figure 9:  Constraints and preferred parameter regions for the CP-violating
Higgs portal model with $\phi=\pi/2$. The dotted lines indicate the ratio of
loop-induced spin-independent interactions and tree-level momentum-suppressed
interactions in terms of the predicted number of events in LZ.

In figure 9 we also indicate the parameter regions excluded by the constraint
$\text{BR}(h\to\text{inv})<0.26$ Aaboud:2018sfi ; Sirunyan:2018owy  as well as
the combinations of $\lambda_{h\chi}/\Lambda$ and $m_{\chi}$ for which the
observed DM relic abundance can be reproduced via annihilations into SM
particles Beniwal:2015sdl . The requirement of EFT validity,
$\lambda_{h\chi}/\Lambda<2\pi/m_{\chi}$ Athron:2018hpc , is satisfied in the
entire parameter region shown in figure 9. We find that for $\phi=\pi/2$
constraints from direct detection experiments are rather weak and only probe
parameter regions where the standard freeze-out calculation predicts $\chi$ to
be a sub-dominant DM component. For these parameter regions we implicitly
assume that the abundance of $\chi$ is set by a non-standard mechanism (e.g. a
particle-antiparticle asymmetry) such that $\chi$ accounts for all of the DM.
If, on the other hand, bounds from direct detection experiments are rescaled
based on the abundance of $\chi$ obtained from standard freeze-out, as done
e.g. in ref. Athron:2018hpc , loop-induced direct detection signals do not
provide relevant constraints on the CP-violating Higgs portal model for the
foreseeable future.

##  4 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 two-loop 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 $\pi/2$ (corresponding to pseudoscalar
interactions). Moreover, they gain in importance as the sensitivity of direct
detection experiments improves. A particularly 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 $\pi/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 $\mathcal{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.

###### Acknowledgements.

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”.

##  Appendix 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 Passarino:1978jh

$\displaystyle\int\frac{\text{d}^{4}k}{(2\pi)^{4}}\frac{1}{[(p+k)^{2}-M^{2}]\,[k^{2}-m^{2}]^{2}}$
$\displaystyle=\frac{i}{(4\pi)^{2}}\,C_{0}(p^{2},\,m^{2},\,M^{2})\;,$ (39)
$\displaystyle\int\frac{\text{d}^{4}k}{(2\pi)^{4}}\frac{k^{\mu}}{[(p+k)^{2}-M^{2}]\,[k^{2}-m^{2}]^{2}}$
$\displaystyle=\frac{i}{(4\pi)^{2}}\,p^{\mu}\,C_{2}(p^{2},\,m^{2},\,M^{2})\;.$
(40)

The $X_{2}$ and $Y_{2}$ functions are given by Abe:2015rja

$\displaystyle\int\frac{\text{d}^{4}k}{(2\pi)^{4}}\frac{1}{[(p+k)^{2}-M^{2}]\,k^{2}\,[k^{2}-m^{2}]^{2}}$
$\displaystyle=\frac{i}{(4\pi)^{2}}\,X_{2}(p^{2},\,M^{2},\,0,\,m^{2})\;,$ (41)
$\displaystyle\int\frac{\text{d}^{4}k}{(2\pi)^{4}}\frac{k^{\mu}}{[(p+k)^{2}-M^{2}]\,k^{2}\,[k^{2}-m^{2}]^{2}}$
$\displaystyle=\frac{i}{(4\pi)^{2}}\,p^{\mu}\,Y_{2}(p^{2},\,M^{2},\,0,\,m^{2})\;,$
(42)

which will reappear in the full two-loop approach in appendix B.1. Finally, we
define the $Z$ functions

$\displaystyle\begin{split}\int\frac{\text{d}^{4}k}{(2\pi)^{4}}\frac{k^{\mu}k^{\nu}}{[(p+k)^{2}-M^{2}]\,k^{4}\,[k^{2}-m^{2}]^{2}}&\\\
&=\frac{i}{(4\pi)^{2}}\Big{(}p^{\mu}p^{\nu}\,Z_{11}(p^{2},\,M^{2},\,m^{2})+g^{\mu\nu}\,Z_{00}(p^{2},\,M^{2},\,m^{2})\Big{)}\;,\end{split}$
(43)
$\displaystyle\begin{split}\int\frac{\text{d}^{4}k}{(2\pi)^{4}}\frac{k^{\mu}k^{\nu}k^{\alpha}}{[(p+k)^{2}-M^{2}]\,k^{4}\,[k^{2}-m^{2}]^{2}}&\\\
=\frac{i}{(4\pi)^{2}}\Big{(}p^{\mu}p^{\nu}p^{\alpha}\,Z_{111}(p^{2},\,M^{2},\,m^{2})&+\big{(}g^{\mu\nu}p^{\alpha}+g^{\alpha\mu}p^{\nu}+g^{\nu\alpha}p^{\mu}\big{)}\,Z_{001}(p^{2},\,M^{2},\,m^{2})\Big{)}\;.\end{split}$
(44)

All of these functions can be calculated readily with Package-X Patel:2015tea
. As examples, we quote the expression for $C_{0}(p^{2},\,m^{2},\,M^{2})$,

$\displaystyle C_{0}(p^{2},\,m^{2},\,M^{2})=$
$\displaystyle-\frac{\log\left(\frac{m^{2}}{M^{2}}\right)}{2p^{2}}+\Bigg{[}\frac{\left(m^{2}-M^{2}-p^{2}\right)}{p^{2}\sqrt{m^{4}-2m^{2}M^{2}-2m^{2}p^{2}+M^{4}-2M^{2}p^{2}+p^{4}}}$
$\displaystyle\times\log\left(\frac{m^{2}+\sqrt{m^{4}-2m^{2}M^{2}-2m^{2}p^{2}+M^{4}-2M^{2}p^{2}+p^{4}}+M^{2}-p^{2}}{2mM}\right)\Bigg{]}\;,$
(45)

and $Z_{11}(M^{2},\,M^{2},\,m^{2})$ evaluated using the on-shell condition
$p^{2}=M^{2}$,

$\displaystyle Z_{11}(M^{2},\,M^{2},\,m^{2})=$
$\displaystyle\,\frac{1}{3m^{2}M^{4}}-\frac{\log\left(\frac{m^{2}}{M^{2}}\right)}{6M^{6}}$
$\displaystyle+\frac{\sqrt{m^{2}\left(m^{2}-4M^{2}\right)}\left(m^{4}-2m^{2}M^{2}-2M^{4}\right)\log\left(\frac{\sqrt{m^{2}\left(m^{2}-4M^{2}\right)}+m^{2}}{2mM}\right)}{3m^{4}M^{6}\left(m^{2}-4M^{2}\right)}\;.$
(46)

The remaining coefficients can be computed analogously.

###  A.2 Box diagram computation and coefficients

Figure 10: Box and corresponding crossed diagram including visualization of
the momentum flow in the limit of no momentum transfer.

For the computation of the box and its crossed diagram shown in figure 10 we
follow the procedure from ref. Abe:2018emu , which allows us to derive the
coefficients for the twist-2 operators. We first start with the amplitude,
which can be expressed as

$\displaystyle
i\mathcal{M}_{\text{Box}}=g_{\chi}^{2}\,g_{\text{SM}}^{2}\,\frac{m_{q}^{2}}{v^{2}}$
$\displaystyle\times\bar{u}_{\chi}(p_{\chi})\int\frac{\text{d}^{4}k}{(2\pi)^{4}}\left(\frac{\not{k}}{(k+p_{\chi})^{2}-m_{\chi}^{2}}+\frac{2\,m_{\chi}\,[\cos^{2}(\phi_{\chi})+i\gamma_{5}\sin(\phi_{\chi})\cos(\phi_{\chi})]}{(k+p_{\chi})^{2}-m_{\chi}^{2}}\right)u_{\chi}(p_{\chi})$
$\displaystyle\times\frac{1}{(k^{2}-m_{a}^{2})^{2}}\,\bar{u}_{q}(p_{q})\left(\frac{-\not{k}}{(k-p_{q})^{2}-m_{q}^{2}}+\frac{2\,m_{q}\,[\cos^{2}(\phi_{\text{SM}})+i\gamma_{5}\sin(\phi_{\text{SM}})\cos(\phi_{\text{SM}})]}{(k-p_{q})^{2}-m_{q}^{2}}\right)u_{q}(p_{q})$
$\displaystyle+\text{crossed diagram}\;,$ (47)

where the crossed diagram is obtained by the replacement $k\rightarrow-k$
within $\bar{u}_{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:

$\displaystyle\frac{-1}{(k-p_{q})^{2}-m_{q}^{2}}+\frac{1}{(k+p_{q})^{2}-m_{q}^{2}}$
$\displaystyle=\frac{-4\,k\cdot p_{q}}{k^{4}}+\mathcal{O}((p_{q})^{2})\;,$
(48)
$\displaystyle\frac{1}{(k-p_{q})^{2}-m_{q}^{2}}+\frac{1}{(k+p_{q})^{2}-m_{q}^{2}}$
$\displaystyle=\frac{2}{k^{2}}+\mathcal{O}((p_{q})^{2})\;,$ (49)

where we have employed the on-shell condition $(p_{q})^{2}=m_{q}^{2}$. The
amplitude then reads

$\displaystyle
i\mathcal{M}_{\text{Box}}=g_{\chi}^{2}\,g_{\text{SM}}^{2}\,\frac{m_{q}^{2}}{v^{2}}$
$\displaystyle\times\bar{u}_{\chi}(p_{\chi})\int\frac{\text{d}^{4}k}{(2\pi)^{4}}\left(\frac{\not{k}}{(k+p_{\chi})^{2}-m_{\chi}^{2}}+\frac{2\,m_{\chi}\,[\cos^{2}(\phi_{\chi})+i\gamma_{5}\sin(\phi_{\chi})\cos(\phi_{\chi})]}{(k+p_{\chi})^{2}-m_{\chi}^{2}}\right)u_{\chi}(p_{\chi})$
$\displaystyle\times\frac{1}{(k^{2}-m_{a}^{2})^{2}}\,\bar{u}_{q}(p_{q})\left(\frac{-\not{k}\,4\,k\cdot
p_{q}}{k^{4}}+\frac{4\,m_{q}\,[\cos^{2}(\phi_{\text{SM}})+i\gamma_{5}\sin(\phi_{\text{SM}})\cos(\phi_{\text{SM}})]}{k^{2}}\right)u_{q}(p_{2})$
$\displaystyle+\mathcal{O}((p_{q})^{2})\;.$ (50)

We can now identify the loop functions defined in appendix A.1, construct the
corresponding effective Lagrangian and use the following decomposition:

$\displaystyle\bar{q}i\partial^{\mu}\gamma^{\nu}q=\mathcal{O}^{q}_{\mu\nu}+\bar{q}\frac{i\partial^{\mu}\gamma^{\nu}-i\partial^{\nu}\gamma^{\mu}}{2}q+\frac{1}{4}g^{\mu\nu}m_{q}\bar{q}q\;.$
(51)

This then yields the effective box diagram Lagrangian given in eqs. (9) and
(10) with the coefficients given by

$\displaystyle\mathcal{C}^{\text{box}}_{1,q}=$
$\displaystyle\,\kappa_{q}\,m_{\chi}\Big{(}-m_{\chi}^{2}\,Z_{111}(m_{\chi}^{2},\,m_{\chi}^{2},\,m_{a}^{2})-6\,Z_{001}(m_{\chi}^{2},\,m_{\chi}^{2},\,m_{a}^{2})$
$\displaystyle+4\,\cos^{2}(\phi_{\text{SM}})\,Y_{2}(m_{\chi}^{2},\,m_{\chi}^{2},\,0,\,m_{a}^{2})-2\,m_{\chi}^{2}\cos^{2}(\phi_{\chi})\,Z_{11}(m_{\chi}^{2},\,m_{\chi}^{2},\,m_{a}^{2})$
(52)
$\displaystyle-8\,\cos^{2}(\phi_{\chi})\,Z_{00}(m_{\chi}^{2},\,m_{\chi}^{2},\,m_{a}^{2})+8\,\cos^{2}(\phi_{\chi})\cos^{2}(\phi_{\text{SM}})\,X_{2}(m_{\chi}^{2},\,m_{\chi}^{2},\,0,\,m_{a}^{2})\Big{)}\;,$
$\displaystyle\mathcal{C}^{\text{box}}_{2,q}=$
$\displaystyle\,\kappa_{q}\,m_{\chi}\Big{(}-2\,m_{\chi}^{2}\cos(\phi_{\chi})\sin(\phi_{\chi})\,Z_{11}(m_{\chi}^{2},\,m_{\chi}^{2},\,m_{a}^{2})$
$\displaystyle-8\,\cos(\phi_{\chi})\sin(\phi_{\chi})\,Z_{00}(m_{\chi}^{2},\,m_{\chi}^{2},\,m_{a}^{2})$
(53)
$\displaystyle+8\,\cos(\phi_{\chi})\sin(\phi_{\chi})\cos^{2}(\phi_{\text{SM}})\,X_{2}(m_{\chi}^{2},\,m_{\chi}^{2},\,0,\,m_{a}^{2})\Big{)}\;,$
$\displaystyle\begin{split}\mathcal{C}^{\text{box}}_{3,q}=&\,\kappa_{q}\,m_{\chi}\Big{(}4\,\cos(\phi_{\text{SM}})\sin(\phi_{\text{SM}})\,Y_{2}(m_{\chi}^{2},\,m_{\chi}^{2},\,0,\,m_{a}^{2})\\\
&+8\,\cos^{2}(\phi_{\chi})\cos(\phi_{\text{SM}})\sin(\phi_{\text{SM}})\,X_{2}(m_{\chi}^{2},\,m_{\chi}^{2},\,0,\,m_{a}^{2})\Big{)}\;,\end{split}$
(54) $\displaystyle\mathcal{C}^{\text{box}}_{4,q}=$
$\displaystyle\,8\,\kappa_{q}\,m_{\chi}\cos(\phi_{\chi})\sin(\phi_{\chi})\cos(\phi_{\text{SM}})\sin(\phi_{\text{SM}})\,X_{2}(m_{\chi}^{2},\,m_{\chi}^{2},\,0,\,m_{a}^{2})\;,$
(55) $\displaystyle\mathcal{C}^{\text{box}}_{5,q}=$
$\displaystyle-8\,\kappa_{q}\,Z_{001}(m_{\chi}^{2},\,m_{\chi}^{2},\,m_{a}^{2})\;,$
(56) $\displaystyle\mathcal{C}^{\text{box}}_{6,q}=$
$\displaystyle\,\kappa_{q}\,m_{\chi}\Big{(}-4\,Z_{111}(m_{\chi}^{2},\,m_{\chi}^{2},\,m_{a}^{2})-8\,\cos^{2}(\phi_{\chi})\,Z_{11}(m_{\chi}^{2},\,m_{\chi}^{2},\,m_{a}^{2})\Big{)}\;,$
(57) $\displaystyle\mathcal{C}^{\text{box}}_{7,q}=$
$\displaystyle-8\,\kappa_{q}\,m_{\chi}\cos(\phi_{\chi})\sin(\phi_{\chi})\,Z_{11}(m_{\chi}^{2},\,m_{\chi}^{2},\,m_{a}^{2})\;,$
(58)

where we have used the shorthand notation

$\displaystyle\kappa_{q}=\frac{g^{2}_{\chi}\,g^{2}_{\text{SM}}\,m_{q}^{2}}{16v^{2}\pi^{2}}\;.$
(59)

##  Appendix 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 Passarino:1978jh

$\displaystyle\int\frac{\text{d}^{4}k}{(2\pi)^{4}}\frac{1}{[(p+k)^{2}-M^{2}]\,[k^{2}-m^{2}]^{3}}$
$\displaystyle=\frac{i}{(4\pi)^{2}}\,D_{0}(p^{2},\,m^{2},\,M^{2})\;,$ (60)
$\displaystyle\int\frac{\text{d}^{4}k}{(2\pi)^{4}}\frac{k^{\mu}}{[(p+k)^{2}-M^{2}]\,[k^{2}-m^{2}]^{3}}$
$\displaystyle=\frac{i}{(4\pi)^{2}}\,p^{\mu}\,D_{3}(p^{2},\,m^{2},\,M^{2})\;.$
(61)

We further define the $X_{n}$ and $Y_{n}$ functions by Abe:2015rja

$\displaystyle\int\frac{\text{d}^{4}k}{(2\pi)^{4}}\frac{1}{[k^{2}-\frac{m_{q}^{2}}{x(1-x)}]^{n}\,[(p_{\chi}+k)^{2}-m_{\chi}^{2}]\,[k^{2}-m_{a}^{2}]}$
$\displaystyle=\frac{i}{(4\pi)^{2}}\,X_{n}\left(p_{\chi}^{2},\,m_{\chi}^{2},\,m_{a}^{2},\,\tfrac{m_{q}^{2}}{x(1-x)}\right)\;,$
(62)
$\displaystyle\int\frac{\text{d}^{4}k}{(2\pi)^{4}}\frac{k^{\mu}}{[k^{2}-\frac{m_{q}^{2}}{x(1-x)}]^{n}\,[(p_{\chi}+k)^{2}-m_{\chi}^{2}]\,[k^{2}-m_{a}^{2}]}$
$\displaystyle=\frac{i}{(4\pi)^{2}}\,p_{\chi}^{\mu}\,Y_{n}\left(p_{\chi}^{2},\,m_{\chi}^{2},\,m_{a}^{2},\,\tfrac{m_{q}^{2}}{x(1-x)}\right)\;.$
(63)

Using partial fraction decomposition, we can derive the following relations
for $X_{n}$

$\displaystyle\begin{split}X_{1}&\left(p_{\chi}^{2},\,m_{\chi}^{2},\,m_{a}^{2},\,\tfrac{m_{q}^{2}}{x(1-x)}\right)\\\
&=\frac{1}{m_{a}^{2}-\frac{m_{q}^{2}}{x(1-x)}}\left(B_{0}\left(p_{\chi}^{2},\,m_{a}^{2},\,m_{\chi}^{2}\right)-B_{0}\left(p_{\chi}^{2},\,\tfrac{m_{q}^{2}}{x(1-x)},\,m_{\chi}^{2}\right)\right)\;,\end{split}$
(64)
$\displaystyle\begin{split}X_{2}&\left(p_{\chi}^{2},\,m_{\chi}^{2},\,m_{a}^{2},\,\tfrac{m_{q}^{2}}{x(1-x)}\right)\\\
&=\frac{1}{m_{a}^{2}-\frac{m_{q}^{2}}{x(1-x)}}\left(X_{1}\left(p_{\chi}^{2},\,m_{\chi}^{2},\,m_{a}^{2},\,\tfrac{m_{q}^{2}}{x(1-x)}\right)-C_{0}\left(p_{\chi}^{2},\,\tfrac{m_{q}^{2}}{x(1-x)},\,m_{\chi}^{2}\right)\right)\;,\end{split}$
(65)
$\displaystyle\begin{split}X_{3}&\left(p_{\chi}^{2},\,m_{\chi}^{2},\,m_{a}^{2},\,\tfrac{m_{q}^{2}}{x(1-x)}\right)\\\
&=\frac{1}{m_{a}^{2}-\frac{m_{q}^{2}}{x(1-x)}}\left(X_{2}\left(p_{\chi}^{2},\,m_{\chi}^{2},\,m_{a}^{2},\,\tfrac{m_{q}^{2}}{x(1-x)}\right)-D_{0}\left(p_{\chi}^{2},\,\tfrac{m_{q}^{2}}{x(1-x)},\,m_{\chi}^{2}\right)\right)\;,\end{split}$
(66)

as well as for $Y_{n}$ Abe:2018emu

$\displaystyle\begin{split}Y_{1}&\left(p_{\chi}^{2},\,m_{\chi}^{2},\,m_{a}^{2},\,\tfrac{m_{q}^{2}}{x(1-x)}\right)\\\
&=\frac{1}{m_{a}^{2}-\frac{m_{q}^{2}}{x(1-x)}}\left(B_{1}\left(p_{\chi}^{2},\,m_{a}^{2},\,m_{\chi}^{2}\right)-B_{1}\left(p_{\chi}^{2},\,\tfrac{m_{q}^{2}}{x(1-x)},\,m_{\chi}^{2}\right)\right)\;,\end{split}$
(67)
$\displaystyle\begin{split}Y_{2}&\left(p_{\chi}^{2},\,m_{\chi}^{2},\,m_{a}^{2},\,\tfrac{m_{q}^{2}}{x(1-x)}\right)\\\
&=\frac{1}{m_{a}^{2}-\frac{m_{q}^{2}}{x(1-x)}}\left(Y_{1}\left(p_{\chi}^{2},\,m_{\chi}^{2},\,m_{a}^{2},\,\tfrac{m_{q}^{2}}{x(1-x)}\right)-C_{2}\left(p_{\chi}^{2},\,\tfrac{m_{q}^{2}}{x(1-x)},\,m_{\chi}^{2}\right)\right)\;,\end{split}$
(68)
$\displaystyle\begin{split}Y_{3}&\left(p_{\chi}^{2},\,m_{\chi}^{2},\,m_{a}^{2},\,\tfrac{m_{q}^{2}}{x(1-x)}\right)\\\
&=\frac{1}{m_{a}^{2}-\frac{m_{q}^{2}}{x(1-x)}}\left(Y_{2}\left(p_{\chi}^{2},\,m_{\chi}^{2},\,m_{a}^{2},\,\tfrac{m_{q}^{2}}{x(1-x)}\right)-D_{3}\left(p_{\chi}^{2},\,\tfrac{m_{q}^{2}}{x(1-x)},\,m_{\chi}^{2}\right)\right)\;,\end{split}$
(69)

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 Hahn:1998yk
.111111We could in principle evaluate the functions $X_{n}$ and $Y_{n}$
directly with Package-X Patel:2015tea  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. We further need the derivatives of these
functions with respect to $m_{a}^{2}$

$\displaystyle\begin{split}\frac{\partial}{\partial
m_{a}^{2}}X_{1}&\left(p_{\chi}^{2},\,m_{\chi}^{2},\,m_{a}^{2},\,\tfrac{m_{q}^{2}}{x(1-x)}\right)\\\
&=\frac{1}{m_{a}^{2}-\frac{m_{q}^{2}}{x(1-x)}}\left(C_{0}(p_{\chi}^{2},\,m_{a}^{2},\,m_{\chi}^{2})-X_{1}\left(p_{\chi}^{2},\,m_{\chi}^{2},\,m_{a}^{2},\,\tfrac{m_{q}^{2}}{x(1-x)}\right)\right)\;,\end{split}$
(70) $\displaystyle\frac{\partial}{\partial m_{a}^{2}}X_{2}$
$\displaystyle\left(p_{\chi}^{2},\,m_{\chi}^{2},\,m_{a}^{2},\,\tfrac{m_{q}^{2}}{x(1-x)}\right)$
$\displaystyle=\frac{1}{m_{a}^{2}-\frac{m_{q}^{2}}{x(1-x)}}\left(\frac{\partial}{\partial
m_{a}^{2}}X_{1}\left(p_{\chi}^{2},\,m_{\chi}^{2},\,m_{a}^{2},\,\tfrac{m_{q}^{2}}{x(1-x)}\right)-X_{2}\left(p_{\chi}^{2},\,m_{\chi}^{2},\,m_{a}^{2},\,\tfrac{m_{q}^{2}}{x(1-x)}\right)\right)\;,$
(71) $\displaystyle\frac{\partial}{\partial m_{a}^{2}}X_{3}$
$\displaystyle\left(p_{\chi}^{2},\,m_{\chi}^{2},\,m_{a}^{2},\,\tfrac{m_{q}^{2}}{x(1-x)}\right)$
$\displaystyle=\frac{1}{m_{a}^{2}-\frac{m_{q}^{2}}{x(1-x)}}\left(\frac{\partial}{\partial
m_{a}^{2}}X_{2}\left(p_{\chi}^{2},\,m_{\chi}^{2},\,m_{a}^{2},\,\tfrac{m_{q}^{2}}{x(1-x)}\right)-X_{3}\left(p_{\chi}^{2},\,m_{\chi}^{2},\,m_{a}^{2},\,\tfrac{m_{q}^{2}}{x(1-x)}\right)\right)\;,$
(72)

as well as Abe:2018emu

$\displaystyle\begin{split}\frac{\partial}{\partial
m_{a}^{2}}Y_{1}&\left(p_{\chi}^{2},\,m_{\chi}^{2},\,m_{a}^{2},\,\tfrac{m_{q}^{2}}{x(1-x)}\right)\\\
&=\frac{1}{m_{a}^{2}-\frac{m_{q}^{2}}{x(1-x)}}\left(C_{2}(p_{\chi}^{2},\,m_{a}^{2},\,m_{\chi}^{2})-Y_{1}\left(p_{\chi}^{2},\,m_{\chi}^{2},\,m_{a}^{2},\,\tfrac{m_{q}^{2}}{x(1-x)}\right)\right)\;,\end{split}$
(73) $\displaystyle\frac{\partial}{\partial m_{a}^{2}}Y_{2}$
$\displaystyle\left(p_{\chi}^{2},\,m_{\chi}^{2},\,m_{a}^{2},\,\tfrac{m_{q}^{2}}{x(1-x)}\right)$
$\displaystyle=\frac{1}{m_{a}^{2}-\frac{m_{q}^{2}}{x(1-x)}}\left(\frac{\partial}{\partial
m_{a}^{2}}Y_{1}\left(p_{\chi}^{2},\,m_{\chi}^{2},\,m_{a}^{2},\,\tfrac{m_{q}^{2}}{x(1-x)}\right)-Y_{2}\left(p_{\chi}^{2},\,m_{\chi}^{2},\,m_{a}^{2},\,\tfrac{m_{q}^{2}}{x(1-x)}\right)\right)\;,$
(74) $\displaystyle\frac{\partial}{\partial m_{a}^{2}}Y_{3}$
$\displaystyle\left(p_{\chi}^{2},\,m_{\chi}^{2},\,m_{a}^{2},\,\tfrac{m_{q}^{2}}{x(1-x)}\right)$
$\displaystyle=\frac{1}{m_{a}^{2}-\frac{m_{q}^{2}}{x(1-x)}}\left(\frac{\partial}{\partial
m_{a}^{2}}Y_{2}\left(p_{\chi}^{2},\,m_{\chi}^{2},\,m_{a}^{2},\,\tfrac{m_{q}^{2}}{x(1-x)}\right)-Y_{3}\left(p_{\chi}^{2},\,m_{\chi}^{2},\,m_{a}^{2},\,\tfrac{m_{q}^{2}}{x(1-x)}\right)\right)\;.$
(75)

###  B.2 Review of the two-loop computation

The two-loop computation was recently presented for a purely pseudoscalar
mediator in the context of a 2HDM Abe:2018emu . 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. Abe:2018emu . 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 Hisano:2010ct ;
Novikov:1983gd .

Figure 11: Visualization of the first step in the two-loop calculation where
we compute the heavy quark loop contribution to the effective vertices between
$a$ and gluons Novikov:1983gd . Here $p^{\prime}$ is given by
$p^{\prime}=p-q$.

The general amplitude can thus be written as

$\displaystyle\begin{split}i\mathcal{M}^{aaG}_{\text{2-Loop}}=-&\sum_{q=c,b,t}\left(\frac{ig_{\text{SM}}\,m_{q}}{v}\right)^{2}\\\
&\times\int\frac{\text{d}^{4}p}{(2\pi)^{4}}\text{Tr}^{\text{Dirac}}_{\text{Colour}}\bigg{[}\big{[}\cos(\phi_{\text{SM}})+i\gamma_{5}\sin(\phi_{\text{SM}})\big{]}\,iS(p)\,\big{[}\cos(\phi_{\text{SM}})+i\gamma_{5}\sin(\phi_{\text{SM}})\big{]}\,i\tilde{S}(p-q)\,\bigg{]}\;,\end{split}$
(76)

where two different coloured fermion propagators, $S(p)$ and $\tilde{S}(p-q)$,
occur, as the Fock-Schwinger gauge breaks translational invariance. Since we
are interested in the effective vertices $aa\,G^{a}_{\mu\nu}G^{a\mu\nu}$ and
$aa\,G^{a}_{\mu\nu}\widetilde{G}^{a\mu\nu}$, we have to consider the following
terms in the propagators121212Our expressions for $S(p)$ and $\tilde{S}(p)$
differ from ref. Hisano:2010ct  by the replacement $G^{a}_{\mu\nu}\rightarrow
G^{a}_{\nu\mu}$. The resulting minus sign drops out, however, as we are only
interested in terms involving two gluon field strength tensors.

$\displaystyle iS(p)=$
$\displaystyle\,iS^{(0)}(p)+\int\text{d}^{4}k_{1}\,iS^{(0)}(p)\,\,g_{s}\gamma^{\alpha}\frac{1}{2}G_{\mu\alpha}\bigg{(}\frac{\partial}{\partial
k_{1\mu}}\delta^{(4)}(k_{1})\bigg{)}\,\,iS^{(0)}(p-k_{1})$
$\displaystyle+\int\text{d}^{4}k_{1}\text{d}^{4}k_{2}\,iS^{(0)}(p)\,\,g_{s}\gamma^{\alpha}\frac{1}{2}G_{\mu\alpha}\bigg{(}\frac{\partial}{\partial
k_{1\mu}}\delta^{(4)}(k_{1})\bigg{)}\,\,iS^{(0)}(p-k_{1})$ (77)
$\displaystyle\hskip 14.22636pt\times
g_{s}\gamma^{\beta}\frac{1}{2}G_{\nu\beta}\bigg{(}\frac{\partial}{\partial
k_{2\nu}}\delta^{(4)}(k_{2})\bigg{)}\,\,iS^{(0)}(p-k_{1}-k_{2})+...\;,$
$\displaystyle i\tilde{S}(p)=$
$\displaystyle\,iS^{(0)}(p)+\int\text{d}^{4}k_{1}\,iS^{(0)}(p+k_{1})\,\,g_{s}\gamma^{\beta}\frac{1}{2}G_{\nu\beta}\bigg{(}\frac{\partial}{\partial
k_{1\nu}}\delta^{(4)}(k_{1})\bigg{)}\,\,iS^{(0)}(p)$
$\displaystyle+\int\text{d}^{4}k_{1}\text{d}^{4}k_{2}\,iS^{(0)}(p+k_{1}+k_{2})\,\,g_{s}\gamma^{\alpha}\frac{1}{2}G_{\mu\alpha}\bigg{(}\frac{\partial}{\partial
k_{2\mu}}\delta^{(4)}(k_{2})\bigg{)}$ (78) $\displaystyle\hskip
14.22636pt\times
iS^{(0)}(p+k_{1})\,\,g_{s}\gamma^{\beta}\frac{1}{2}G_{\nu\beta}\bigg{(}\frac{\partial}{\partial
k_{1\nu}}\delta^{(4)}(k_{1})\bigg{)}\,\,iS^{(0)}(p)+...\;,$

where terms with derivatives acting on $G_{\mu\nu}$ are neglected since they
are not relevant for the present work. Further we have used
$iS^{(0)}(p)=i(\not{p}+m)/(p^{2}-m^{2})$ and $G^{\mu\nu}=G^{a\mu\nu}\,t^{a}$
with $t^{a}$ being an $SU(3)$ generator fulfilling
$\text{Tr}[t^{a}t^{b}]=\delta^{ab}/2$. Inserting these expressions into eq.
(76), we identify three terms relevant for the computation of vertices
involving two gluon field strength tensors. These terms are visualised in
figure 11. We will present the various calculational steps for the term
involving a gluon field strength tensor each from $S(p)$ and $\tilde{S}(p-q)$,
illustrated in the left panel in figure 11. The corresponding term reads

$\displaystyle i\mathcal{M}^{aaG}_{\text{2-Loop}}\supset-$
$\displaystyle\sum_{q=c,b,t}\left(\frac{ig_{\text{SM}}\,m_{q}}{v}\right)^{2}$
$\displaystyle\times\frac{g_{s}^{2}}{8}\,G^{a}_{\alpha\mu}G^{a}_{\beta\nu}\int\frac{\text{d}^{4}p}{(2\pi)^{4}}\frac{\partial}{\partial
k_{1\mu}}\frac{\partial}{\partial
k_{2\nu}}\text{Tr}^{\text{Dirac}}\bigg{[}\big{[}\cos(\phi_{\text{SM}})+i\gamma_{5}\sin(\phi_{\text{SM}})\big{]}\frac{\not{p}+m}{p^{2}-m^{2}}\gamma^{\alpha}\frac{\not{p}-\not{k}_{1}+m}{(p-k_{1})^{2}-m^{2}}$
$\displaystyle\times\big{[}\cos(\phi_{\text{SM}})+i\gamma_{5}\sin(\phi_{\text{SM}})\big{]}\frac{\not{p}-\not{q}+\not{k}_{2}+m}{(p-q+k_{2})^{2}-m^{2}}\gamma^{\beta}\frac{\not{p}-\not{q}+m}{(p-q)^{2}-m^{2}}\bigg{]}_{k_{1}=k_{2}=0}\;,$
(79)

where we performed the trace over colour indices as well as partial
integration regarding $k_{1}$ and $k_{2}$ and switched the indices of the
gluon field strengths. After performing the derivatives with respect to
$k_{1\mu}$ and $k_{2\nu}$, setting $k_{1}=k_{2}=0$ afterwards and evaluating
the trace using Package-X Patel:2015tea , we now have to project out the
leading spin-independent and spin-dependent term. The spin-independent term
can be obtained by rewriting131313The signs in front of
$\mathcal{O}^{G}_{\mu\nu}$ and $\mathcal{O}^{G}_{\alpha\beta}$ differ from eq.
(50) in ref. Hisano:2010ct . As we are not including the terms containing the
twist-2 operator, this difference will not play a role.

$\displaystyle\begin{split}G^{a}_{\alpha\mu}G^{a}_{\beta\nu}&=\frac{1}{12}G^{a}_{\rho\sigma}G^{a\rho\sigma}(g_{\alpha\beta}g_{\mu\nu}-g_{\alpha\nu}g_{\beta\mu})+\frac{1}{2}g_{\alpha\beta}\,\mathcal{O}^{G}_{\mu\nu}\\\
&\quad+\frac{1}{2}g_{\mu\nu}\,\mathcal{O}^{G}_{\alpha\beta}-\frac{1}{2}g_{\alpha\nu}\,\mathcal{O}^{G}_{\beta\mu}-\frac{1}{2}g_{\beta\mu}\,\mathcal{O}^{G}_{\alpha\nu}+\mathcal{O}^{G}_{\alpha\mu\beta\nu}\;,\end{split}$
(80)

where we have introduced the twist-2 gluon operator $\mathcal{O}^{G}_{\mu\nu}$
and a higher spin operator $\mathcal{O}^{G}_{\alpha\mu\beta\nu}$

$\displaystyle\mathcal{O}^{G}_{\mu\nu}$
$\displaystyle=G^{a\rho}_{\hphantom{a\rho}\mu}G^{a}_{\hphantom{a}\rho\nu}-\frac{1}{4}g_{\mu\nu}G^{a}_{\rho\sigma}G^{a\rho\sigma}\;,$
(81)
$\displaystyle\begin{split}\mathcal{O}^{G}_{\alpha\mu\beta\nu}&=G^{a}_{\alpha\mu}G^{a}_{\beta\nu}-\frac{1}{2}g_{\alpha\beta}G^{a\rho}_{\hphantom{a\rho}\mu}G^{a}_{\hphantom{a}\rho\nu}-\frac{1}{2}g_{\mu\nu}G^{a\rho}_{\hphantom{a\rho}\alpha}G^{a}_{\hphantom{a}\rho\beta}\\\
&\quad+\frac{1}{2}g_{\alpha\nu}G^{a\rho}_{\hphantom{a\rho}\beta}G^{a}_{\hphantom{a}\rho\mu}+\frac{1}{2}g_{\beta\mu}G^{a\rho}_{\hphantom{a\rho}\alpha}G^{a}_{\hphantom{a}\rho\nu}\\\
&\quad+\frac{1}{6}G^{a}_{\rho\sigma}G^{a\rho\sigma}(g_{\alpha\beta}g_{\mu\nu}-g_{\alpha\nu}g_{\beta\mu})\;.\end{split}$
(82)

These operators do not contribute to $G^{a}_{\rho\sigma}G^{a\rho\sigma}$ and
give sub-leading SI interactions such that they are neglected in the present
work. We therefore only have to consider the replacement

$\displaystyle
G^{a}_{\alpha\mu}G^{a}_{\beta\nu}\rightarrow\frac{1}{12}G^{a}_{\rho\sigma}G^{a\rho\sigma}(g_{\alpha\beta}g_{\mu\nu}-g_{\alpha\nu}g_{\beta\mu})\;,$
(83)

for the spin-independent term. The spin-dependent term involving
$G^{a}_{\mu\nu}\widetilde{G}^{a\mu\nu}$ can be obtained straightforwardly by
isolating the term in the trace containing the $\epsilon$-tensor. Putting both
contributions together, we then obtain

$\displaystyle
i\mathcal{M}^{aaG}_{\text{2-Loop}}\supset-\sum_{q=c,b,t}\bigg{(}$
$\displaystyle\frac{ig_{\text{SM}}\,m_{q}}{v}\bigg{)}^{2}\Bigg{\\{}\frac{g_{s}^{2}}{96}\,G^{a}_{\rho\sigma}G^{a\rho\sigma}\int\frac{\text{d}^{4}p}{(2\pi)^{4}}\,\frac{-48\,[2\,m_{q}^{2}\,\cos(2\phi_{\text{SM}})+p\cdot(p-q)]}{[p^{2}-m_{q}^{2}]^{2}\,[(p-q)^{2}-m_{q}^{2}]^{2}}$
$\displaystyle+g_{s}^{2}\,G^{a}_{\rho\sigma}\widetilde{G}^{a\rho\sigma}\int\frac{\text{d}^{4}p}{(2\pi)^{4}}\frac{m_{q}^{2}\,\sin(2\phi_{\text{SM}})}{[p^{2}-m_{q}^{2}]^{2}\,[(p-q)^{2}-m_{q}^{2}]^{2}}\Bigg{\\}}\;.$
(84)

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:

$\displaystyle\frac{1}{A^{2}\,B^{2}}=\int_{0}^{1}\text{d}x\,\frac{6\,x(1-x)}{(xA+(1-x)B)^{4}}\;,$
(85) $\displaystyle
x((p-q)^{2}-m^{2})+(1-x)(p^{2}-m^{2})=(p-qx)^{2}-m^{2}-q^{2}x(-1+x)\;.$ (86)

Shifting $p\rightarrow p+qx$ and performing the loop integral over $p$, we
obtain finally

$\displaystyle i\mathcal{M}^{aaG}_{\text{2-Loop}}\supset$
$\displaystyle\,i\sum_{q=c,b,t}\left(\frac{g_{\text{SM}}\,m_{q}}{v}\right)^{2}\frac{g^{2}_{s}}{32\pi^{2}}\,G^{a}_{\rho\sigma}G^{a\rho\sigma}$
$\displaystyle\hskip
35.56593pt\times\int_{0}^{1}\text{d}x\,\Bigg{\\{}\frac{-2\,m_{q}^{2}\,[\cos(2\phi_{\text{SM}})-\frac{1}{2}]\,x(1-x)}{[m_{q}^{2}-q^{2}x(1-x)]^{2}}+\frac{x(1-x)}{[m_{q}^{2}-q^{2}x(1-x)]}\Bigg{\\}}$
$\displaystyle\hskip
2.84544pt+i\sum_{q=c,b,t}\left(\frac{g_{\text{SM}}\,m_{q}}{v}\right)^{2}\frac{g^{2}_{s}}{16\pi^{2}}\,G^{a}_{\rho\sigma}{\widetilde{G}}^{a\rho\sigma}\int_{0}^{1}\text{d}x\,\frac{m_{q}^{2}\,\sin(2\phi_{\text{SM}})\,x(1-x)}{[m_{q}^{2}-q^{2}x(1-x)]^{2}}\;.$
(87)

One can proceed in a similar fashion for the remaining two terms in eq. (76),
which contribute equally to the amplitude, such that the full amplitude reads

$\displaystyle i\mathcal{M}^{aaG}_{\text{2-Loop}}=$
$\displaystyle\,i\,d^{\text{\,full}}_{G}(q^{2})\frac{\alpha_{s}}{12\pi}\,G^{a}_{\rho\sigma}G^{a\rho\sigma}+i\,d^{\text{\,full}}_{\widetilde{G}}(q^{2})\frac{\alpha_{s}}{8\pi}\,G^{a}_{\rho\sigma}\widetilde{G}^{a\rho\sigma}\;,$
(88)

with

$\displaystyle\begin{split}d^{\text{\,full}}_{G}(q^{2})=&\,\sum_{q=c,b,t}\left(\frac{g_{\text{SM}}\,m_{q}}{v}\right)^{2}\int_{0}^{1}\text{d}x\,\Bigg{\\{}\frac{\frac{3}{2}\,x(1-x)}{[m_{q}^{2}-q^{2}x(1-x)]}\\\
&+\frac{m_{q}^{2}}{2}\,\frac{3\,(1-x)x+2\,(-1-x+x^{2})\,\cos(2\phi_{\text{SM}})}{[m_{q}^{2}-q^{2}x(1-x)]^{2}}-m_{q}^{4}\,\frac{1-3x+3x^{2}-(1-x)x\,\cos(2\phi_{\text{SM}})}{[m_{q}^{2}-q^{2}x(1-x)]^{3}}\Bigg{\\}}\;,\end{split}$
(89) $\displaystyle d^{\text{\,full}}_{\widetilde{G}}(q^{2})=$
$\displaystyle\,\sum_{q=c,b,t}\left(\frac{g_{\text{SM}}\,m_{q}}{v}\right)^{2}\int_{0}^{1}\text{d}x\,\frac{m_{q}^{2}\,\sin(2\phi_{\text{SM}})}{[m_{q}^{2}-q^{2}x(1-x)]^{2}}\;,$
(90)

where we have made the integral over the Feynman parameter symmetric under
$x\rightarrow 1-x$. This result agrees with eq. (B.73) from ref. Abe:2018emu
for $\phi_{\text{SM}}=\pi/2$. Moreover, when performing a heavy quark
expansion on eq. (88), we recover the effective Lagrangian from eq. (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

$\displaystyle\begin{split}i\mathcal{M}^{\text{SI}}_{\text{2-Loop}}=-\bar{u}_{\chi}(p_{\chi})\,g_{\chi}^{2}\int&\frac{\text{d}^{4}q}{(2\pi)^{4}}\,\frac{\not{q}+2\,m_{\chi}\,[\cos^{2}(\phi_{\chi})+i\gamma_{5}\sin(\phi_{\chi})\cos(\phi_{\chi})]}{[(p_{\chi}+q)^{2}-m_{\chi}^{2}]\,[q^{2}-m_{a}^{2}]^{2}}\,u_{\chi}(p_{\chi})\\\
&\times\bigg{[}d^{\text{\,full}}_{G}(q^{2})\,\frac{\alpha_{s}}{12\pi}\,G^{a}_{\rho\sigma}G^{a\rho\sigma}+d^{\text{\,full}}_{\widetilde{G}}(q^{2})\,\frac{\alpha_{s}}{8\pi}\,G^{a}_{\rho\sigma}\widetilde{G}^{a\rho\sigma}\bigg{]}\;.\end{split}$
(91)

In terms of the various loop functions defined in appendix B.1, we can map
this amplitude onto

$\displaystyle\mathcal{L}^{\text{SI}}_{\text{2-Loop}}$
$\displaystyle=\left(\mathcal{C}^{\text{full}}_{G,S}\,\bar{\chi}\chi+\mathcal{C}^{\text{full}}_{G,PS}\,\bar{\chi}i\gamma_{5}\chi\,\right)\frac{-\alpha_{s}}{12\pi}\,G^{a}_{\mu\nu}G^{a\mu\nu}\;,$
(92) $\displaystyle\mathcal{L}^{\text{SD}}_{\text{2-Loop}}$
$\displaystyle=\left(\mathcal{C}^{\text{full}}_{\widetilde{G},S}\,\bar{\chi}\chi+\mathcal{C}^{\text{full}}_{\widetilde{G},PS}\,\bar{\chi}i\gamma_{5}\chi\,\right)\frac{\alpha_{s}}{8\pi}\,G^{a}_{\mu\nu}{\widetilde{G}}^{a\mu\nu}\;,$
(93)

with

$\displaystyle\begin{split}\mathcal{C}^{\text{full}}_{G,S}=\frac{1}{(4\pi)^{2}}\sum_{q=c,b,t}&\left(\frac{g_{\chi}\,g_{\text{SM}}\,m_{q}}{v}\right)^{2}\\\
&\times\bigg{\\{}m_{\chi}\,F_{Y}(p_{\chi}^{2},\,m_{\chi}^{2},\,m_{a}^{2},\,m_{q}^{2})+2\,m_{\chi}\cos^{2}(\phi_{\chi})\,F_{X}(p_{\chi}^{2},\,m_{\chi}^{2},\,m_{a}^{2},\,m_{q}^{2})\bigg{\\}}\;,\end{split}$
(94)
$\displaystyle\mathcal{C}^{\text{full}}_{G,PS}=\frac{1}{(4\pi)^{2}}\sum_{q=c,b,t}$
$\displaystyle\left(\frac{g_{\chi}\,g_{\text{SM}}\,m_{q}}{v}\right)^{2}2\,m_{\chi}\sin(\phi_{\chi})\cos(\phi_{\chi})\,F_{X}(p_{\chi}^{2},\,m_{\chi}^{2},\,m_{a}^{2},\,m_{q}^{2})\;,$
(95)
$\displaystyle\begin{split}\mathcal{C}^{\text{full}}_{\widetilde{G},S}=-\frac{1}{(4\pi)^{2}}\sum_{q=c,b,t}&\left(\frac{g_{\chi}\,g_{\text{SM}}\,m_{q}}{v}\right)^{2}\\\
&\times\bigg{\\{}m_{\chi}\,E_{Y}(p_{\chi}^{2},\,m_{\chi}^{2},\,m_{a}^{2},\,m_{q}^{2})+2\,m_{\chi}\cos^{2}(\phi_{\chi})\,E_{X}(p_{\chi}^{2},\,m_{\chi}^{2},\,m_{a}^{2},\,m_{q}^{2})\bigg{\\}}\;,\end{split}$
(96)
$\displaystyle\mathcal{C}^{\text{full}}_{\widetilde{G},PS}=-\frac{1}{(4\pi)^{2}}\sum_{q=c,b,t}$
$\displaystyle\left(\frac{g_{\chi}\,g_{\text{SM}}\,m_{q}}{v}\right)^{2}2\,m_{\chi}\sin(\phi_{\chi})\cos(\phi_{\chi})\,E_{X}(p_{\chi}^{2},\,m_{\chi}^{2},\,m_{a}^{2},\,m_{q}^{2})\;.$
(97)

Here we have introduced the shorthand notation

$\displaystyle
F_{\Lambda}(p_{\chi}^{2},\,m_{\chi}^{2},\,m_{a}^{2},\,m_{q}^{2})$
$\displaystyle=\int_{0}^{1}\text{d}x\bigg{[}-\frac{3}{2}\,\frac{\partial}{\partial
m_{a}^{2}}\,\Lambda_{1}\left(p_{\chi}^{2},\,m_{\chi}^{2},\,m_{a}^{2},\,\tfrac{m_{q}^{2}}{x(1-x)}\right)$
$\displaystyle+\frac{m_{q}^{2}}{2}\,\frac{3\,(1-x)x+2\,(-1-x+x^{2})\,\cos(2\phi_{\text{SM}})}{x^{2}(1-x)^{2}}\frac{\partial}{\partial
m_{a}^{2}}\,\Lambda_{2}\left(p_{\chi}^{2},\,m_{\chi}^{2},\,m_{a}^{2},\,\tfrac{m_{q}^{2}}{x(1-x)}\right)$
$\displaystyle+m_{q}^{4}\,\frac{1-3x+3x^{2}-(1-x)x\,\cos(2\phi_{\text{SM}})}{x^{3}(1-x)^{3}}\frac{\partial}{\partial
m_{a}^{2}}\,\Lambda_{3}\left(p_{\chi}^{2},\,m_{\chi}^{2},\,m_{a}^{2},\,\tfrac{m_{q}^{2}}{x(1-x)}\right)\bigg{]}\;,$
(98) $\displaystyle
E_{\Lambda}(p_{\chi}^{2},\,m_{\chi}^{2},\,m_{a}^{2},\,m_{q}^{2})$
$\displaystyle=\int_{0}^{1}\text{d}x\,\bigg{[}m_{q}^{2}\,\frac{\sin(2\phi_{\text{SM}})}{x^{2}(1-x)^{2}}\,\frac{\partial}{\partial
m_{a}^{2}}\,\Lambda_{2}\left(p_{\chi}^{2},\,m_{\chi}^{2},\,m_{a}^{2},\,\tfrac{m_{q}^{2}}{x(1-x)}\right)\bigg{]}\;,$
(99)

with $\Lambda=X,Y$.

##  Appendix 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 Shifman:1978zn ;
Jungman:1995df :

$\displaystyle q$ $\displaystyle=u,d,s:$ $\displaystyle\quad\langle
N|m_{q}\bar{q}q|N\rangle$ $\displaystyle=m_{N}f^{N}_{q}\;,$ (100)
$\displaystyle Q$ $\displaystyle=c,b,t:$ $\displaystyle\quad\langle
N|m_{Q}\bar{Q}Q|N\rangle$ $\displaystyle=\langle
N|-\frac{\alpha_{s}}{12\pi}G^{a}_{\mu\nu}G^{a\mu\nu}|N\rangle=\frac{2}{27}m_{N}f^{N}_{G}\;,$
(101) $\displaystyle q$ $\displaystyle=u,d,s,c,b:$ $\displaystyle\quad\langle
N|\mathcal{O}^{q}_{\mu\nu}|N\rangle$
$\displaystyle=\frac{1}{m_{N}}\Big{(}p^{N}_{\mu}p^{N}_{\nu}-\frac{1}{4}m_{N}^{2}g_{\mu\nu}\Big{)}\Big{(}q^{N}(2)+\bar{q}^{N}(2)\Big{)}\;,$
(102)

where $f^{N}_{q}$ and $f^{N}_{G}$ are form factors, $m_{N}$ is the nucleon
mass, $q^{N}(2)$ and $\bar{q}^{N}(2)$ are the second moments of the quark
parton distribution functions and $p^{N}_{\mu}$ is the nucleon four-momentum.
The values of the form factors for light quarks are taken from micrOmegas
Belanger:2018mqt 141414We refer to ref. Ellis:2018dmb  for a discussion of the
uncertainties of these form factors, in particular regarding the strange quark
matrix element.

$\displaystyle f^{p}_{u}$ $\displaystyle=0.01513\;,\qquad$ $\displaystyle
f^{p}_{d}$ $\displaystyle=0.0191\;,\qquad$ $\displaystyle f^{p}_{s}$
$\displaystyle=0.0447\;,$ (103) $\displaystyle f^{n}_{u}$
$\displaystyle=0.0110\;,\qquad$ $\displaystyle f^{n}_{d}$
$\displaystyle=0.0273\;,\qquad$ $\displaystyle f^{n}_{s}$
$\displaystyle=0.0447\;,$ (104)

which are related to the gluon form factor via Jungman:1995df

$\displaystyle f^{p}_{G}=1-\sum_{q=u,d,s}f^{p}_{q}=0.92107\;,\qquad
f^{n}_{G}=1-\sum_{q=u,d,s}f^{n}_{q}=0.917\;.$ (105)

The second moments are calculated at the scale $\mu=m_{Z}$ by using CTEQ PDFs
Abe:2018emu ; Pumplin:2002vw . For the proton, one finds

$\displaystyle u^{p}(2)$ $\displaystyle=0.22\;,\qquad$
$\displaystyle\bar{u}^{p}(2)$ $\displaystyle=0.034\;,$ (106) $\displaystyle
d^{p}(2)$ $\displaystyle=0.11\;,\qquad$ $\displaystyle\bar{d}^{p}(2)$
$\displaystyle=0.036\;,$ (107) $\displaystyle s^{p}(2)$
$\displaystyle=0.026\;,\qquad$ $\displaystyle\bar{s}^{p}(2)$
$\displaystyle=0.026\;,$ (108) $\displaystyle c^{p}(2)$
$\displaystyle=0.019\;,\qquad$ $\displaystyle\bar{c}^{p}(2)$
$\displaystyle=0.019\;,$ (109) $\displaystyle b^{p}(2)$
$\displaystyle=0.012\;,\qquad$ $\displaystyle\bar{b}^{p}(2)$
$\displaystyle=0.012\;,$ (110)

whereas for the neutron the up- and down-quark values have to be interchanged.

For the spin-dependent interactions we need the following form factors:

$\displaystyle q$ $\displaystyle=u,d,s:$ $\displaystyle\quad\langle
N^{\prime}|m_{q}\bar{q}i\gamma_{5}q|N\rangle$
$\displaystyle=F^{q/N}_{P}(q^{2})\;,$ (111) $\displaystyle Q$
$\displaystyle=c,b,t:$ $\displaystyle\quad\langle
N^{\prime}|m_{Q}\bar{Q}i\gamma_{5}Q|N\rangle$ $\displaystyle=\langle
N^{\prime}|\frac{\alpha_{s}}{8\pi}G^{a}_{\mu\nu}{\widetilde{G}}^{a\mu\nu}|N\rangle=F^{N}_{\widetilde{G}}(q^{2})\;,$
(112)

where $N^{\prime}$ refers to a change of nucleon momentum. This explicit
dependence on the momentum transfer $q^{\mu}$ arises from non-negligible $\pi$
and $\eta$ pole contributions. The corresponding expressions are given in eqs.
(A30) and (A42) of ref. Bishara:2017pfq  and are implemented in DirectDM
Bishara:2017nnn .

## References

  * (1) P. Cushman et al., Working Group Report: WIMP Dark Matter Direct Detection,  in Proceedings, 2013 Community Summer Study on the Future of U.S. Particle Physics, 2013.  1310.8327.
  * (2) PandaX-II, J. Xia et al., Constraining WIMP-Nucleon Effective Interactions from PandaX-II Experiment, 1807.01936.
  * (3) A. L. Fitzpatrick, W. Haxton, E. Katz, N. Lubbers, and Y. Xu, The Effective Field Theory of Dark Matter Direct Detection,  JCAP 1302 (2013) 004, [1203.3542].
  * (4) N. Anand, A. L. Fitzpatrick, and W. C. Haxton, Weakly interacting massive particle-nucleus elastic scattering response,  Phys. Rev. C89 (2014), no. 6 065501, [1308.6288].
  * (5) R. Catena and P. Gondolo, Global fits of the dark matter-nucleon effective interactions,  JCAP 1409 (2014), no. 09 045, [1405.2637].
  * (6) M. I. Gresham and K. M. Zurek, Effect of nuclear response functions in dark matter direct detection,  Phys. Rev. D89 (2014), no. 12 123521, [1401.3739].
  * (7) R. Catena, Prospects for direct detection of dark matter in an effective theory approach,  JCAP 1407 (2014) 055, [1406.0524].
  * (8) V. Gluscevic, M. I. Gresham, S. D. McDermott, A. H. G. Peter, and K. M. Zurek, Identifying the Theory of Dark Matter with Direct Detection,  JCAP 1512 (2015), no. 12 057, [1506.04454].
  * (9) J. B. Dent, L. M. Krauss, J. L. Newstead, and S. Sabharwal, General analysis of direct dark matter detection: From microphysics to observational signatures,  Phys. Rev. D92 (2015), no. 6 063515, [1505.03117].
  * (10) F. Kahlhoefer and S. Wild, Studying generalised dark matter interactions with extended halo-independent methods,  JCAP 1610 (2016), no. 10 032, [1607.04418].
  * (11) F. Bishara, J. Brod, B. Grinstein, and J. Zupan, Chiral Effective Theory of Dark Matter Direct Detection,  JCAP 1702 (2017), no. 02 009, [1611.00368].
  * (12) T. D. P. Edwards, B. J. Kavanagh, and C. Weniger, Assessing Near-Future Direct Dark Matter Searches with Benchmark-Free Forecasting,  Phys. Rev. Lett. 121 (2018), no. 18 181101, [1805.04117].
  * (13) U. Haisch and F. Kahlhoefer, On the importance of loop-induced spin-independent interactions for dark matter direct detection,  JCAP 1304 (2013) 050, [1302.4454].
  * (14) A. Crivellin, F. D’Eramo, and M. Procura, New Constraints on Dark Matter Effective Theories from Standard Model Loops,  Phys. Rev. Lett. 112 (2014) 191304, [1402.1173].
  * (15) A. Crivellin and U. Haisch, Dark matter direct detection constraints from gauge bosons loops,  Phys. Rev. D90 (2014) 115011, [1408.5046].
  * (16) F. D’Eramo, B. J. Kavanagh, and P. Panci, You can hide but you have to run: direct detection with vector mediators,  JHEP 08 (2016) 111, [1605.04917].
  * (17) F. Bishara, J. Brod, B. Grinstein, and J. Zupan, Renormalization Group Effects in Dark Matter Interactions, 1809.03506.
  * (18) M. Freytsis and Z. Ligeti, On dark matter models with uniquely spin-dependent detection possibilities,  Phys. Rev. D83 (2011) 115009, [1012.5317].
  * (19) K. R. Dienes, J. Kumar, B. Thomas, and D. Yaylali, Overcoming Velocity Suppression in Dark-Matter Direct-Detection Experiments,  Phys. Rev. D90 (2014), no. 1 015012, [1312.7772].
  * (20) C. Boehm, M. J. Dolan, C. McCabe, M. Spannowsky, and C. J. Wallace, Extended gamma-ray emission from Coy Dark Matter,  JCAP 1405 (2014) 009, [1401.6458].
  * (21) S. Ipek, D. McKeen, and A. E. Nelson, A Renormalizable Model for the Galactic Center Gamma Ray Excess from Dark Matter Annihilation,  Phys. Rev. D90 (2014), no. 5 055021, [1404.3716].
  * (22) J. M. No, Looking through the pseudoscalar portal into dark matter: Novel mono-Higgs and mono-Z signatures at the LHC,  Phys. Rev. D93 (2016), no. 3 031701, [1509.01110].
  * (23) D. Goncalves, P. A. N. Machado, and J. M. No, Simplified Models for Dark Matter Face their Consistent Completions,  Phys. Rev. D95 (2017), no. 5 055027, [1611.04593].
  * (24) M. Bauer, U. Haisch, and F. Kahlhoefer, Simplified dark matter models with two Higgs doublets: I. Pseudoscalar mediators,  JHEP 05 (2017) 138, [1701.07427].
  * (25) M. Bauer, M. Klassen, and V. Tenorth, Universal properties of pseudoscalar mediators in dark matter extensions of 2HDMs,  JHEP 07 (2018) 107, [1712.06597].
  * (26) P. Pani and G. Polesello, Dark matter production in association with a single top-quark at the LHC in a two-Higgs-doublet model with a pseudoscalar mediator,  Phys. Dark Univ. 21 (2018) 8–15, [1712.03874].
  * (27) P. Tunney, J. M. No, and M. Fairbairn, Probing the pseudoscalar portal to dark matter via $\bar{b}bZ(\to\ell\ell)+\not{E}_{T}$ : From the LHC to the Galactic Center excess,  Phys. Rev. D96 (2017), no. 9 095020, [1705.09670].
  * (28) S. Banerjee, D. Barducci, G. Bélanger, B. Fuks, A. Goudelis, et al., Cornering pseudoscalar-mediated dark matter with the LHC and cosmology,  JHEP 07 (2017) 080, [1705.02327].
  * (29) U. Haisch, J. F. Kamenik, A. Malinauskas, and M. Spira, Collider constraints on light pseudoscalars,  JHEP 03 (2018) 178, [1802.02156].
  * (30) LHC Dark Matter Working Group, T. Abe et al., LHC Dark Matter Working Group: Next-generation spin-0 dark matter models, 1810.09420.
  * (31) B. Batell, M. Pospelov, and A. Ritz, Exploring Portals to a Hidden Sector Through Fixed Targets,  Phys. Rev. D80 (2009) 095024, [0906.5614].
  * (32) M. Freytsis, Z. Ligeti, and J. Thaler, Constraining the Axion Portal with $B\to Kl^{+}l^{-}$,  Phys. Rev. D81 (2010) 034001, [0911.5355].
  * (33) B. Batell, M. Pospelov, and A. Ritz, Multi-lepton Signatures of a Hidden Sector in Rare B Decays,  Phys. Rev. D83 (2011) 054005, [0911.4938].
  * (34) M. J. Dolan, F. Kahlhoefer, C. McCabe, and K. Schmidt-Hoberg, A taste of dark matter: Flavour constraints on pseudoscalar mediators,  JHEP 03 (2015) 171, [1412.5174]. [Erratum: JHEP07,103(2015)].
  * (35) A. Berlin, S. Gori, T. Lin, and L.-T. Wang, Pseudoscalar Portal Dark Matter,  Phys. Rev. D92 (2015) 015005, [1502.06000].
  * (36) B. Döbrich, F. Ertas, F. Kahlhoefer, and T. Spadaro, Model-independent bounds on light pseudoscalars from rare B-meson decays,  Phys. Lett. B790 (2019) 537–544, [1810.11336].
  * (37) J. Billard, L. Strigari, and E. Figueroa-Feliciano, Implication of neutrino backgrounds on the reach of next generation dark matter direct detection experiments,  Phys. Rev. D89 (2014), no. 2 023524, [1307.5458].
  * (38) G. Arcadi, M. Lindner, F. S. Queiroz, W. Rodejohann, and S. Vogl, Pseudoscalar Mediators: A WIMP model at the Neutrino Floor,  JCAP 1803 (2018), no. 03 042, [1711.02110].
  * (39) N. F. Bell, G. Busoni, and I. W. Sanderson, Loop Effects in Direct Detection,  JCAP 1808 (2018), no. 08 017, [1803.01574].
  * (40) T. Li, Revisiting the direct detection of dark matter in simplified models,  Phys. Lett. B782 (2018) 497–502, [1804.02120].
  * (41) T. Abe, M. Fujiwara, and J. Hisano, Loop corrections to dark matter direct detection in a pseudoscalar mediator dark matter model, 1810.01039.
  * (42) LUX-ZEPLIN, D. S. Akerib et al., Projected WIMP Sensitivity of the LUX-ZEPLIN (LZ) Dark Matter Experiment, 1802.06039.
  * (43) XENON, E. Aprile et al., Physics reach of the XENON1T dark matter experiment,  JCAP 1604 (2016), no. 04 027, [1512.07501].
  * (44) F. Kahlhoefer, K. Schmidt-Hoberg, and S. Wild, Dark matter self-interactions from a general spin-0 mediator,  JCAP 1708 (2017), no. 08 003, [1704.02149].
  * (45) A. Beniwal, F. Rajec, C. Savage, P. Scott, C. Weniger, et al., Combined analysis of effective Higgs portal dark matter models,  Phys. Rev. D93 (2016), no. 11 115016, [1512.06458].
  * (46) GAMBIT, P. Athron et al., Global analyses of Higgs portal singlet dark matter models using GAMBIT, 1808.10465.
  * (47) T. Abe and R. Sato, Current status and future prospects of the singlet-doublet dark matter model with CP-violation, 1901.02278.
  * (48) G. D’Ambrosio, G. F. Giudice, G. Isidori, and A. Strumia, Minimal flavor violation: An Effective field theory approach,  Nucl. Phys. B645 (2002) 155–187, [hep-ph/0207036].
  * (49) J. Abdallah et al., Simplified Models for Dark Matter Searches at the LHC,  Phys. Dark Univ. 9-10 (2015) 8–23, [1506.03116].
  * (50) J. Ellis, TikZ-Feynman: Feynman diagrams with TikZ,  Comput. Phys. Commun. 210 (2017) 103–123, [1601.05437].
  * (51) J. Hisano, K. Ishiwata, and N. Nagata, Gluon contribution to the dark matter direct detection,  Phys. Rev. D82 (2010) 115007, [1007.2601].
  * (52) J. Hisano, R. Nagai, and N. Nagata, Effective Theories for Dark Matter Nucleon Scattering,  JHEP 05 (2015) 037, [1502.02244].
  * (53) M. A. Shifman, A. I. Vainshtein, and V. I. Zakharov, Remarks on Higgs Boson Interactions with Nucleons,  Phys. Lett. 78B (1978) 443–446.
  * (54) H. H. Patel, Package-X: A Mathematica package for the analytic calculation of one-loop integrals,  Comput. Phys. Commun. 197 (2015) 276–290, [1503.01469].
  * (55) J. Fan, M. Reece, and L.-T. Wang, Non-relativistic effective theory of dark matter direct detection,  JCAP 1011 (2010) 042, [1008.1591].
  * (56) The GAMBIT Dark Matter Workgroup, T. Bringmann et al., DarkBit: A GAMBIT module for computing dark matter observables and likelihoods,  Eur. Phys. J. C77 (2017), no. 12 831, [1705.07920].
  * (57) F. Bishara, J. Brod, B. Grinstein, and J. Zupan, DirectDM: a tool for dark matter direct detection, 1708.02678.
  * (58) XENON, E. Aprile et al., Dark Matter Search Results from a One Ton-Year Exposure of XENON1T,  Phys. Rev. Lett. 121 (2018), no. 11 111302, [1805.12562].
  * (59) PandaX-II, X. Cui et al., Dark Matter Results From 54-Ton-Day Exposure of PandaX-II Experiment,  Phys. Rev. Lett. 119 (2017), no. 18 181302, [1708.06917].
  * (60) LUX, D. S. Akerib et al., Results from a search for dark matter in the complete LUX exposure,  Phys. Rev. Lett. 118 (2017), no. 2 021303, [1608.07648].
  * (61) LUX, D. S. Akerib et al., Limits on spin-dependent WIMP-nucleon cross section obtained from the complete LUX exposure,  Phys. Rev. Lett. 118 (2017), no. 25 251302, [1705.03380].
  * (62) M. Kaplinghat, S. Tulin, and H.-B. Yu, Direct Detection Portals for Self-interacting Dark Matter,  Phys. Rev. D89 (2014), no. 3 035009, [1310.7945].
  * (63) G. Belanger, F. Boudjema, A. Goudelis, A. Pukhov, and B. Zaldivar, micrOMEGAs5.0 : Freeze-in,  Comput. Phys. Commun. 231 (2018) 173–186, [1801.03509].
  * (64) BaBar, J. P. Lees et al., Search for di-muon decays of a low-mass Higgs boson in radiative decays of the Υ(1S),  Phys. Rev. D87 (2013), no. 3 031102, [1210.0287]. [Erratum: Phys. Rev.D87,no.5,059903(2013)].
  * (65) BaBar, J. P. Lees et al., Search for hadronic decays of a light Higgs boson in the radiative decay $\Upsilon\to\gamma A^{0}$,  Phys. Rev. Lett. 107 (2011) 221803, [1108.3549].
  * (66) W. Altmannshofer, P. Paradisi, and D. M. Straub, Model-Independent Constraints on New Physics in $b\to s$ Transitions,  JHEP 04 (2012) 008, [1111.1257].
  * (67) LHCb, R. Aaij et al., Measurement of the $B^{0}_{s}\to\mu^{+}\mu^{-}$ branching fraction and effective lifetime and search for $B^{0}\to\mu^{+}\mu^{-}$ decays,  Phys. Rev. Lett. 118 (2017), no. 19 191801, [1703.05747].
  * (68) C. Bobeth, M. Gorbahn, T. Hermann, M. Misiak, E. Stamou, et al., $B_{s,d}\to l^{+}l^{-}$ in the Standard Model with Reduced Theoretical Uncertainty,  Phys. Rev. Lett. 112 (2014) 101801, [1311.0903].
  * (69) C.-Y. Chen, H. Davoudiasl, W. J. Marciano, and C. Zhang, Implications of a light “dark Higgs” solution to the $g_{μ}$-2 discrepancy,  Phys. Rev. D93 (2016), no. 3 035006, [1511.04715].
  * (70) W. J. Marciano, A. Masiero, P. Paradisi, and M. Passera, Contributions of axionlike particles to lepton dipole moments,  Phys. Rev. D94 (2016), no. 11 115033, [1607.01022].
  * (71) S. Mantry, M. Pitschmann, and M. J. Ramsey-Musolf, Distinguishing axions from generic light scalars using electric dipole moment and fifth-force experiments,  Phys. Rev. D90 (2014), no. 5 054016, [1401.7339].
  * (72) ATLAS, CMS, G. Aad et al., Measurements of the Higgs boson production and decay rates and constraints on its couplings from a combined ATLAS and CMS analysis of the LHC pp collision data at $\sqrt{s}=7$ and 8 TeV,  JHEP 08 (2016) 045, [1606.02266].
  * (73) ATLAS, M. Aaboud et al., Search for invisible Higgs boson decays in vector boson fusion at $\sqrt{s}=13$ TeV with the ATLAS detector, Submitted to: Phys. Lett. (2018) [1809.06682].
  * (74) CMS, A. M. Sirunyan et al., Search for invisible decays of a Higgs boson produced through vector boson fusion in proton-proton collisions at $\sqrt{s}=$ 13 TeV, Submitted to: Phys. Lett. (2018) [1809.05937].
  * (75) G. Passarino and M. J. G. Veltman, One Loop Corrections for e+ e- Annihilation Into mu+ mu- in the Weinberg Model,  Nucl. Phys. B160 (1979) 151–207.
  * (76) T. Abe and R. Sato, Quantum corrections to the spin-independent cross section of the inert doublet dark matter,  JHEP 03 (2015) 109, [1501.04161].
  * (77) T. Hahn and M. Perez-Victoria, Automatized one loop calculations in four-dimensions and D-dimensions,  Comput. Phys. Commun. 118 (1999) 153–165, [hep-ph/9807565].
  * (78) V. A. Novikov, M. A. Shifman, A. I. Vainshtein, and V. I. Zakharov, Calculations in External Fields in Quantum Chromodynamics. Technical Review, Fortsch. Phys. 32 (1984) 585.
  * (79) G. Jungman, M. Kamionkowski, and K. Griest, Supersymmetric dark matter,  Phys. Rept. 267 (1996) 195–373, [hep-ph/9506380].
  * (80) J. Ellis, N. Nagata, and K. A. Olive, Uncertainties in WIMP Dark Matter Scattering Revisited,  Eur. Phys. J. C78 (2018), no. 7 569, [1805.09795].
  * (81) J. Pumplin, D. R. Stump, J. Huston, H. L. Lai, P. M. Nadolsky, et al., New generation of parton distributions with uncertainties from global QCD analysis,  JHEP 07 (2002) 012, [hep-ph/0201195].
  * (82) F. Bishara, J. Brod, B. Grinstein, and J. Zupan, From quarks to nucleons in dark matter direct detection,  JHEP 11 (2017) 059, [1707.06998].