# Neutron star observations of pseudoscalar-mediated dark matter

John Coffey<EMAIL_ADDRESS>Department of Physics and Astronomy, University
of Victoria, Victoria, BC V8P 5C2, Canada  David McKeen<EMAIL_ADDRESS>TRIUMF, 4004 Wesbrook Mall, Vancouver, BC V6T 2A3, Canada  David E. Morrissey
<EMAIL_ADDRESS>TRIUMF, 4004 Wesbrook Mall, Vancouver, BC V6T 2A3, Canada
Nirmal Raj<EMAIL_ADDRESS>TRIUMF, 4004 Wesbrook Mall, Vancouver, BC V6T 2A3,
Canada

###### Abstract

Scattering interactions between dark matter and Standard Model states mediated
by pseudoscalars are generically challenging to uncover at direct detection
experiments due to rates suppressed by powers of the local dark matter
velocity $v_{\rm DM}\sim 10^{-3}\,c$. However, they may be observed in the
dark matter-induced heating of neutron stars, whose steep gravitational
potentials prevent such suppression by accelerating infalling particles to
semi-relativistic speeds. We investigate this phenomenon in the context of two
specific, self-consistent scenarios for pseudoscalars coupled to dark matter,
and compare the sensitivity of neutron star heating to bounds from direct
searches for the mediators and dark matter. The first “lighter” scenario
consists of sub-10 GeV mass dark matter mediated by an axion-like particle
(ALP), while the second “heavier” scenario has dark matter above
$10\,\text{GeV}$ mediated by a dark pseudoscalar that mixes with a
pseudoscalar from a two-Higgs doublet (the so-called 2HDM+$a$ model). In both
frameworks, we show that imminent measurements of neutron stars will be able
to test pseudoscalar-mediated dark matter beyond the reach of direct dark
matter searches as well as bounds on the mediators from flavor observables,
beam dump experiments, and high-energy colliders.

##  I Introduction

Figure 1:  Sensitivities to the ALP-gluon coupling versus dark matter mass, of
neutron star heating from capture of DM mediated by the ALP (solid) and the
radial mode $r$ (dashed) in our benchmark scenario, for ALP-to-DM mass ratios
of 1/10 (left) and 1 (right). The bands represent uncertainties in the mass-
radius configurations of candidate neutron stars. Also shown are limits from
dark matter direct searches (brown region) and those on the mediator from beam
dumps, rare meson decays, and astrophysics (gray region). For further details,
see Sec. V.

Figure 2: Sensitivities to the angle of mixing between the pseudoscalars in
the 2HDM+$a$ scenario versus dark matter mass, of neutron star heating from
capture of DM mediated by the pseudoscalar $a$ for our chosen benchmark point,
for $a$-to-DM mass ratios of 1/3 (left) and 1 (right). The bands represent
uncertainties in the mass-radius configurations of candidate neutron stars.
Also shown are limits from dark matter direct detection (brown region),
contours of cross sections within the reach of future detectors, and limits on
the 2HDM+$a$ framework from colliders and flavor observables (gray region).
For further details, see Sec. V.

Even as the hunt for dark matter (DM) is continuing on an unprecedented scale,
its identity remains a closely guarded secret. Direct searches for DM in the
form of a weakly interacting massive particle (WIMP) Goodman and Witten
(1985); Jungman  _et al._ (1996) have placed remarkable limits on DM over the
mass range of $m_{\text{DM}}\sim 1\,\text{GeV}$–$100\,\text{TeV}$ Akerib  _et
al._ (2017); Aprile  _et al._ (2018); Meng  _et al._ (2021). The absence of
WIMP signals in direct detection experiments has also motivated searches for
other DM candidates over a broader range of masses Mitridate  _et al._ (2022);
Adhikari  _et al._ (2022); Carney  _et al._ (2022). Notably, direct detection
collaborations have made substantial progress in the “light” DM regime of sub-
GeV masses, where the lower target recoil energies make scattering more
difficult to detect, with current limits down to $m_{\text{DM}}\sim
10\,\text{MeV}$ in nuclear recoils and $m_{\text{DM}}\sim 100\,\text{keV}$ in
electron recoils Kahn and Lin (2022).

Dark matter could also have avoided direct detection in other ways. For
instance, should scattering of DM on targets proceed inelastically to an
excited state (presumably present in the theory due to an approximate
symmetry), the target recoil rates diminish rapidly with increasing mass
splitting between the dark states Tucker-Smith and Weiner (2001, 2005).
Likewise, the scattering cross section of local DM in a target can be highly
suppressed if the underlying dynamics make it proportional to positive powers
of the DM velocity, $v_{\text{DM}}\sim 10^{-3}$ Fan  _et al._ (2010); Freytsis
and Ligeti (2011); Fitzpatrick  _et al._ (2013); Dienes  _et al._ (2014). This
latter scenario is realized when DM connects to visible matter through a
pseudoscalar mediator. A comprehensive program to discover DM via scattering
must necessarily include probes beyond direct detection.

In recent times, signals from DM capture in compact stars have emerged as a
highly promising avenue; for reviews, see Refs. Tinyakov  _et al._ (2021);
Berti  _et al._ (2022). In particular, capture in neutron stars (NS),
beginning with Ref. Goldman and Nussinov (1989), has attracted much attention
by virtue of properties favorable to DM capture, namely their steep
gravitational potentials and high densities. One of us recently showed with
other collaborators that old, isolated NSs would heat appreciably by the
transfer of kinetic energy from infalling DM, and that the phenomenon could be
detected with forthcoming infrared telescopes trained at nearby NSs Baryakhtar
_et al._ (2017). The observation of a population of anomalously warm NSs would
therefore provide evidence for DM, while the discovery of even a single
relatively cold NS could be used to exclude a broad range of DM candidates.
Much theoretical work has followed this proposal Bramante  _et al._ (2017);
Raj  _et al._ (2018); Acevedo  _et al._ (2020); Chen and Lin (2018); Bell  _et
al._ (2018a); Garani  _et al._ (2019); Camargo  _et al._ (2019); Hamaguchi
_et al._ (2019a); Keung  _et al._ (2020); Bell  _et al._ (2020); Dasgupta  _et
al._ (2020); Garani  _et al._ (2021); Maity and Queiroz (2021); Bell  _et al._
(2021a); Zeng  _et al._ (2022); Anzuini  _et al._ (2021); Bell  _et al._
(2019); Garani and Heeck (2019); Joglekar  _et al._ (2020a, b); Bell  _et al._
(2021b); Fujiwara  _et al._ (2022); Hamaguchi  _et al._ (2022), exploring
several important particle and nuclear astrophysics aspects.

NS heating by dark matter is also strongly complementary to other search
methods. In particular, this indirect probe is sensitive to many types of
particle-like DM that are challenging to detect through direct scattering in
the laboratory. This includes lighter sub-GeV DM Bell  _et al._ (2019);
Joglekar  _et al._ (2020a, b); Acevedo  _et al._ (2020); Keung  _et al._
(2020), DM that scatters primarily inelastically Baryakhtar  _et al._ (2017);
Bell  _et al._ (2018a); Acevedo  _et al._ (2020); Fujiwara  _et al._ (2022),
and scenarios where the scattering rate of DM with nuclei is suppressed by
positive powers of the DM velocity Raj  _et al._ (2018); Bell  _et al._
(2019); Joglekar  _et al._ (2020b).

In this study, we examine the sensitivity of the NS kinetic heating mechanism
to DM that connects to the SM through a pseudoscalar mediator. As mentioned
above, pseudoscalar mediation suppresses scattering cross sections in direct
detection by multiple powers of DM velocity. Neutron star heating overcomes
this suppression by virtue of the semi-relativistic speeds with which DM falls
into NSs. This fact has been noted before (in, e.g. Refs. Raj  _et al._
(2018); Bell  _et al._ (2019); Joglekar  _et al._ (2020b)), but only in the
context of a generic pseudoscalar in the “heavy mediator” limit, i.e. in an
effective field theory (EFT) regime where the mediator mass exceeds the
momentum transfer. We expand on this previous work by investigating the
sensitivity of NSs to pseudoscalar-mediated DM within self-consistent,
ultraviolet (UV) complete scenarios. This allows us to extend beyond the range
of EFT validity to cases where the mediator is light relative to the relevant
4-momentum exchange, and to compare the reach of NS heating observations to
other tests of the DM particle as well as the mediator itself.

To be concrete, in this work we study fermionic dark matter $\chi$ and two
realizations of lighter pseudoscalar mediators. The first realization is in
the form of an axion-like particle (ALP) that couples primarily to gluons, for
which a UV completion may be specified, and on which constraints have been
placed by a multitude of experiments. Since the ALP can be naturally light
relative to the weak scale, we focus on ALP mediators in the context of
lighter DM with $m_{\chi}\lesssim 10\,\text{GeV}$. The second mediator
realization we consider is the well-studied Two Higgs Doublet Model plus
pseudoscalar, the 2HDM+$a$. Here, the lighter pseudoscalar relevant to DM
detection couples primarily to quarks and has a mass that is typically near
the weak scale. We focus on this scenario to study DM in the traditional WIMP
mass range, $m_{\chi}\gtrsim 10$ GeV.

In both these scenarios we find regions where the discovery potential of NS
heating is a considerable improvement over existing constraints. Our primary
results are summarized in Fig. 1 for the ALP mediator and Fig. 2 for the
2HDM+$a$ mediator. ALPs are generically excluded by beam dump searches for
masses below a few 100 MeV, but NS heating would probe surviving islands of
parameter space in this region as well as higher masses. The 2HDM+$a$ scenario
for certain parametric ranges remarkably survives a host of collider and
flavor constraints in the so-called alignment limit, and we identify the
regions that NS heating can exclusively probe.

This paper is organized as follows. Following this introduction, we describe
our benchmark ALP and 2HDM+$a$ scenarios in Section II. Next, in Section III
we provide expressions for nucleon scattering cross sections in both
scenarios. In Section IV we review neutron star capture of dark matter with an
emphasis on kinetic heating, and discuss how these apply specifically to
pseudoscalar mediation. In Section V we present our results. Section VI is
reserved for conclusions and outlook. Various technical details underlying our
main results are presented in the Appendices A, B, C.

##  II Benchmark Scenarios

Our investigation focuses on Dirac fermion dark matter $\chi$ that connects to
the SM primarily through a lighter pseudoscalar mediator. We do not attempt to
specify the cosmological origin of the $\chi$ abundance, but we note that it
can be consistent with bounds from indirect detection if its abundance is
asymmetric Nussinov (1985); Barr  _et al._ (1990); Kaplan  _et al._ (2009).
Furthermore, our assumption of a pseudoscalar mediator also prevents large
non-perturbative enhancements of the annihilation and self-scattering cross
sections, even when the mediator is parametrically lighter than the DM state
Agrawal  _et al._ (2020).

In our study we consider two benchmark scenarios for the pseudoscalar
mediator. In the first, we investigate light (i.e. mass well below the weak
scale) mediators that can arise as axion-like particles (ALPs) from the
spontaneous breaking of an approximate global symmetry. In the second
benchmark, we use the 2HDM+$a$ model in which a singlet pseudoscalar couples
to dark matter directly and connects to the SM by mixing with a heavier
pseudoscalar from a two-Higgs doublet (2HDM) extension of the SM. This second
scenario motivates weak-scale mediator masses. Thus, taking $m_{a}\leq
m_{\chi}$, we focus on the ALP-mediator benchmark for lighter DM with
$m_{\chi}\lesssim 10\,\text{GeV}$ and the 2HDM+$a$ -mediator benchmark for
heavier DM with $m_{\chi}\gtrsim 10\,\text{GeV}$.

###  II.1 Lighter Mediator: Axion-Like Particles

A natural way to realize a light pseudoscalar below the weak scale is in the
form of an axion-like particle (ALP) that emerges as the pseudo-Nambu-
Goldstone boson from the spontaneous breaking of an approximate but anomalous
global $U(1)$ symmetry Peccei and Quinn (1977a, b); Weinberg (1978); Wilczek
(1978). Its mass is set by either the soft breaking scale of the $U(1)$ or the
mass scale of the gauge sector containing the anomaly, and can be
parametrically small. Our first benchmark for a pseudoscalar mediator for DM
consists of an ALP that couples primarily to gluons.

####  II.1.1 Model and Elementary Couplings

In this context we consider a light pseudoscalar arising as the phase of a
complex scalar $\phi$ that couples to DM according to

$\displaystyle\mathcal{L}\ \supset\
-g_{\chi}\,\phi\,\bar{\chi}_{R}\chi_{L}+{\rm h.c.}$ (1)

This is consistent with a global $U(1)_{PQ}$ symmetry under which
$[\phi]_{PQ}=1$ and $[\chi_{L}]_{PQ}-[\chi_{R}]_{PQ}=-1$. If $\phi$ develops a
vacuum expection value (VEV), we can write

$\displaystyle\phi=\frac{f}{\sqrt{2}}\,(1+r/f)\,e^{ia/f}\,$ (2)

with the real scalar $r(x)$ and pseudoscalar ALP $a(x)$. This leads to a mass
for $\chi$,

$\displaystyle m_{\chi}=g_{\chi}\,f/\sqrt{2}\ .$ (3)

Changing variables such that $\chi_{L}\to e^{-i\beta_{L}\,a/f}\chi_{L}$ and
$\chi_{R}\to e^{-i\beta_{R}\,a/f}\chi_{R}$ with $(\beta_{L}-\beta_{R})=1$
removes the pseudoscalar from Eq. (1) at the cost of introducing new
interactions Di Luzio  _et al._ (2020):

$\displaystyle\mathcal{L}$ $\displaystyle\supset$
$\displaystyle\frac{\partial_{\mu}a}{f}\,\big{[}\beta_{L}\,\bar{\chi}_{L}\gamma^{\mu}\chi_{L}+\beta_{R}\,\bar{\chi}_{R}\gamma^{\mu}\chi_{R}\big{]}$
$\displaystyle=$
$\displaystyle\frac{\partial_{\mu}a}{2f}\,\big{[}(\beta_{L}+\beta_{R})\,\bar{\chi}\gamma^{\mu}\chi-\bar{\chi}\gamma^{\mu}\gamma^{5}\chi\big{]}$
$\displaystyle=$
$\displaystyle(\text{E.O.M.})+{m_{\chi}}\,\frac{a}{f}\;\bar{\chi}i\gamma^{5}\chi\
,$

where (E.O.M.) refers to terms that vanish on-shell when the equations of
motion are satisfied. We also include an explicit mass term for the ALP $a(x)$
that softly breaks the global $U(1)$ and allows us to treat the ALP mass
$m_{a}$ as a free parameter.

The pseudoscalar ALP can connect to the SM in the presence of fermions with SM
gauge charges that also get their mass from the VEV of $\phi$. To be concrete,
we assume a set of $N_{Q}$ fermions $Q$ with SM charges $({\bf 3,1,0})$
together with the coupling Kim (1979); Shifman  _et al._ (1980)

$\displaystyle\mathcal{L}\supset-y_{Q}\,\phi\,\bar{Q}_{R}{Q}_{L}+{\rm h.c.}$
(5)

As before, the interaction generates fermion masses $m_{Q}=y_{Q}\,f/\sqrt{2}$
as well as couplings of the form of Eq. (II.1.1) after shifting field
variables. The change of field variables also generates an interaction between
the pseudoscalar and gluons Bauer  _et al._ (2017a, 2021a):

$\displaystyle\mathcal{L}$ $\displaystyle\supset$
$\displaystyle\frac{\alpha_{s}}{4\pi}\,\frac{N_{Q}}{2}\,\frac{a}{f}\,G_{\mu\nu}^{a}\widetilde{G}^{a\,\mu\nu}$
$\displaystyle\equiv$ $\displaystyle
c_{GG}\,\frac{\alpha_{s}}{4\pi}\,\frac{a}{f}\,G_{\mu\nu}^{a}\widetilde{G}^{a\,\mu\nu}\
,$

with $c_{GG}=N_{Q}/2$. We note that the ALP literature often expresses this
coupling in terms of $f_{a}=f/2c_{GG}$ Di Luzio  _et al._ (2020). The
interaction of Eq. (II.1.1) is the primary connection between the ALP and the
SM.

Additional couplings to the SM will be generated by running from the matching
scale $f$ down to the energies relevant for DM scattering. Below the
electroweak scale, these take the form Bauer  _et al._ (2017a, 2021a, 2021b,
2021c)

$\displaystyle\mathcal{L}\supset\frac{\partial_{\mu}a}{f}\big{(}\bar{f}_{L}\gamma^{\mu}\mathbf{k}_{F}f_{L}+\bar{f}_{R}\gamma^{\mu}\mathbf{k_{f}}f_{R}\big{)}$
(7)

where $f_{L,R}$ refer to SM fermions and $\mathbf{k_{F}}$ and $\mathbf{k}_{f}$
are matrices in flavor space. For the dark matter analysis to follow, we will
be interested mainly in the flavor-diagonal coupling combinations
$c_{ii}=(k_{f_{ii}}-k_{F_{ii}})$. Relative to $c_{GG}$ defined in Eq.
(II.1.1), we find $c_{ii}\sim 10^{-2}c_{GG}$.

Together with the ALP, this model also contains the radial mode $r(x)$ of the
complex scalar, $\phi(x)=(f/\sqrt{2})(1+r/f)e^{ia/f}$. In contrast to the ALP,
however, this state obtains a mass directly from spontaneous symmetry breaking
with $m_{r}\propto f$. The most important couplings of $r$ for our analysis
are

$\displaystyle\mathcal{L}\ $ $\displaystyle\supset$ $\displaystyle\
\frac{r}{f}\,(\partial a)^{2}-m_{\chi}\,\frac{r}{f}\,\bar{\chi}\,\chi$
$\displaystyle~{}~{}~{}+\frac{\alpha_{s}}{6\pi}\,c_{GG}\,\frac{r}{f}\,G_{\mu\nu}^{a}G^{a\,\mu\nu}\
,$

where the last term is generated by integrating out the heavy color-charged
fermions for $m_{Q}\gg m_{r}$.

####  II.1.2 Limits on the ALP Model

The most stringent bounds on the theory typically come from the ALP mediator
itself. Experimental limits on ALP-type pseudoscalars coupling to the SM
primarily through the operator of Eq. (II.1.1) have been studied recently in
Refs. Bauer  _et al._ (2017a, 2021a, 2021b, 2021c), while astrophysical limits
from supernova physics are calculated in Refs. Ertas and Kahlhoefer (2020);
Caputo  _et al._ (2022a, b). We review and summarize these limits in App. A
and show the exclusions as shaded regions in our sensitivity plots.

In addition to the ALP, the theory also contains the heavier radial mode $r$
and mediator quarks $Q$ with masses near $m\sim f$. These are discussed in
App. A and can lead to further bounds on the theory, particularly for ALP
masses above $m_{a}\gtrsim 2\,\text{GeV}$. While the detailed bounds depend on
the specific UV details of the model, we find that all current limits can be
satisfied over the entire ALP mass region of interest for $c_{GG}/f\gtrsim
10^{-2}\,\text{GeV}^{-1}$.

###  II.2 Heavier Mediator: 2HDM+$a$

As a benchmark model of a heavier pseudoscalar mediator ($m_{a}\gtrsim
10\,\text{GeV}$), we consider the well-studied 2HDM+$a$ theory consisting of a
two-Higgs doublet model (2HDM) together with an additional pseudoscalar
$a_{0}$ that couples to fermionic dark matter $\chi$ Ipek  _et al._ (2014).

####  II.2.1 Model and Elementary Couplings

In this model a gauge-singlet pseudoscalar $a_{0}$ interacts with the Dirac
fermion DM field $\chi$ via

$\displaystyle\mathcal{L}_{\chi}=y_{\chi}\,a_{0}\,\bar{\chi}i\gamma^{5}\chi\
.$ (9)

The $a_{0}$ also connects to the SM by coupling to the $H_{1}$ and $H_{2}$
doublets of the 2HDM sector according to

$\displaystyle\mathcal{L}_{a12}=-i\,\mu_{a12}\,a_{0}\,H_{1}^{\dagger}H_{2}+{\rm
h.c.}$ (10)

Together, the interactions of Eqs. (9,10) connect the DM particle to the SM
via the pseudoscalar $a_{0}$.

The scalar sector of the 2HDM+$a$ model is described by the potential

$\displaystyle
V=V_{\text{2HDM}}+\frac{1}{2}m_{a_{0}}^{2}a_{0}^{2}-\mathcal{L}_{a12}\ ,$ (11)

where $V_{\text{2HDM}}$ is a standard two-Higgs doublet potential Gunion  _et
al._ (2000); Djouadi (2008). On its own, the 2HDM gives rise to physical
scalars $h$ and $H$, a pseudoscalar $A_{0}$, and a charged complex scalar
$H^{\pm}$. The 2HDM parameters in this potential can be defined in terms of
the expectation values $v_{1}$ and $v_{2}$ of the doublets $H_{1}$ and $H_{2}$
with $v=\sqrt{v_{1}^{2}+v_{2}^{2}}=246\,\text{GeV}$, $\tan\beta=v_{2}/v_{1}$,
the CP-even scalar mixing angle $\alpha$, the 2HDM masses $m_{h}$, $m_{H}$,
$m_{A_{0}}$, $m_{H^{\pm}}$, and the quartic coupling $\lambda_{3}$ connecting
$|H_{1}|^{2}$ and $|H_{2}|^{2}$. We assume further that the $H_{1}$ and
$H_{2}$ doublets couple to the SM fermions in the Type-II format where $H_{1}$
connects only to $d$\- and $e$-type fermions and $H_{2}$ couples exclusively
to $u$-type fermions. Including the term of Eq. (10) leads to mixing between
the $a_{0}$ and $A_{0}$ states described by the angle

$\displaystyle\tan 2\theta=-\frac{2\mu_{a12}v}{m_{A}^{2}-m_{a}^{2}}\ ,$ (12)

where $m_{a}$ and $m_{A}$ are the physical masses of the lighter and heavier
pseudoscalars. Together with the angle $\theta$, these three parameters can be
used to specify the theory in place of $m_{a_{0}}$, $m_{A_{0}}$, and
$\mu_{a12}$.

After electroweak symmetry breaking, the physical pseudoscalar mass
eigenstates $a$ and $A$ couple to DM and SM fermions according to

$\displaystyle-\mathcal{L}_{\rm Yuk}$ $\displaystyle\supset$
$\displaystyle\sum_{f}\frac{m_{f}}{v}(\xi_{a}^{f}a+\xi_{A}^{f}A)\bar{f}i\gamma^{5}f$
(13)
$\displaystyle~{}~{}+(\xi_{a}^{\chi}a+\xi_{A}^{\chi}A)\bar{\chi}i\gamma^{5}\chi$

where for a Type-II 2HDM,

$\displaystyle\begin{array}[]{rclcrcl}\xi_{a}^{u}&=&\sin\theta\,\cot\beta\
,&{}\hfil&\xi_{a}^{d}&=&\sin\theta\,\tan\beta\ ,\\\
\xi_{A}^{u}&=&-\cos\theta\,\cot\beta\
,&{}\hfil&\xi_{A}^{d}&=&-\cos\theta\,\tan\beta\ ,\\\
\xi_{a}^{\chi}&=&y_{\chi}\cos\theta\
,&{}\hfil&\xi_{A}^{\chi}&=&y_{\chi}\sin\theta\ .\end{array}$ (17)

####  II.2.2 Limits on the 2HDM+$a$ Model

For the dark matter studies to follow, we will focus on the model in the limit
of $m_{a}\ll m_{A}$ and $|\theta|\ll 1$ in which the light pseudoscalar is
mostly singlet. The bounds on the theory in this limit then mostly factorize
into those on the 2HDM sector and those on the light pseudoscalar.

As a fixed benchmark for the 2HDM part of the theory, we set $\tan\beta=3$,
$m_{H}=m_{H^{\pm}}=m_{A}=600\,\text{GeV}$, $\lambda_{3}\simeq
m_{h}^{2}/v^{2}$, and work in the alignment/decoupling limit of
$\cos(\beta-\alpha)=0$ such that the lighter $h$ scalar couples to the rest of
the SM in the same way as the SM Higgs boson Gunion  _et al._ (2000); Djouadi
(2008). As discussed further in Appendix B, in the absence of pseudoscalar
mixing these parameters are consistent with direct and indirect bounds on the
Type-II 2HDM including flavour, precision electroweak and Higgs measurements,
and collider searches. It also satisfies the requirements of vacuum stability
and perturbative unitarity Kling  _et al._ (2016).

Within this 2HDM-sector benchmark we allow for variations in $m_{a}$ and
$\theta$. We briefly summarize the bounds on this parameter subspace here, and
discuss them in more detail in Appendix B.

For $m_{a}\lesssim 10\,\text{GeV}$ flavor constraints arise. For
$1\,\text{GeV}\lesssim m_{a}\lesssim 10\,\text{GeV}$ the most important of
these come from the loop-induced $b$-$s$-$a$ coupling Arcadi  _et al._ (2018);
Arnan  _et al._ (2017). This effective coupling leads to rare decays of bottom
mesons, such as $B_{s}\to\mu^{+}\mu^{-}$ and
$B^{\pm}\to{K^{\pm}}^{(\ast)}\mu^{+}\mu^{-}$. Limits from the latter process
Aaij  _et al._ (2015) constrain the mixing angle to be smaller than ${\cal
O}({\rm few}\times 10^{-3})$ for $m_{a}\lesssim 4\,\text{GeV}$.

When $m_{a}\leq m_{h}/2$ the decay channel $h\to aa$ can have a significant
branching fraction for $|\theta|\gtrsim 0.05$ leading to bounds on $\theta$ of
this order from searches for exotic Higgs boson decays as well as Higgs rate
measurements. The strongest limits for larger $m_{a}\geq 90\,\text{GeV}$ come
from LHC searches for $H/A\to\tau\tau$ and imply $|\theta|\lesssim 0.1$ up to
$m_{a}\simeq 2\,m_{t}$ (for $\tan\beta=3$). These bounds are shown in Fig. 6
in App. B and are summarized by the shaded grey regions in Fig. 2.

##  III Dark Matter Direct Detection

(a) spin- and velocity-dependent via exchange of ALP & pion, or 2HDM+$a$
pseudoscalars.

(b) spin- and velocity-independent via exchange of ALP radial mode.

(c) spin- and velocity-independent loops involving 2HDM+$a$ pseudoscalars.

Figure 3: Feynman diagrams leading to dark matter scattering on nucleons in
our setups. The mediators involved in our ALP and 2HDM+$a$ pseudoscalar
scenarios are respectively labeled in red and blue.

In this section we provide expressions for the non-relativistic nuclear
scattering cross sections for DM direct detection in the two benchmark models
presented above. Both models give rise to spin-dependent scattering on nuclei
through direct pseudoscalar exchange as well as subdominant spin-independent
scattering by other mechanisms. We present the leading contributions to each
channel.

###  III.1 ALP Model Results

The dominant contributions to spin-dependent and spin-independent scattering
in the ALP scenario are shown in panels (a) and (b) of Fig. 3, respectively.
The former arises from exchanging the pseudoscalar $a$ together with its
mixing with the neutral pion, and the latter from the radial mode $r(x)$ in
the $t$-channel. Both are dominated by the couplings of the scalars to gluons
rather than to quarks.

####  III.1.1 ALP Spin-Dependent

The interactions of Eqs. (II.1.1,7) connect the axion to nucleons and lead to
mixing with the neutral mesons. Thus, the axion serves as a mediator between
DM and nucleons. To estimate these effects, we work in the limit of two quark
flavors as in Refs. Bauer  _et al._ (2021b, c) (although see Ref. Blinov  _et
al._ (2022) for a related analysis that includes three flavors). We estimate
that an explicit three-flavor treatment leads to only small changes at low
momentum exchange $Q\equiv\sqrt{|q^{2}|}\lesssim 500\,\text{MeV}$, while for
higher momenta other uncertainties are expected to be more important.

The mixing between the axion and the neutral pion in this context depends on
the small parameter

$\displaystyle\varepsilon\equiv\frac{f_{\pi}}{2f}\left[c_{uu}-c_{dd}+2c_{GG}\left(\frac{m_{d}-m_{u}}{m_{u}+m_{d}}\right)\right]\
,$ (18)

where $f_{\pi}\simeq 93\,\text{MeV}$ is the pion decay constant. To leading
non-trivial order in $\varepsilon$, axion-pion mixing is described by the
mixing angle $\alpha$ given by

$\displaystyle\tan(2\alpha)=-2\varepsilon\,\frac{\text{min}\\{m_{a}^{2},\,m_{\pi}^{2}\\}}{m_{\pi}^{2}-m_{a}^{2}}\
,$ (19)

with physical masses

$\displaystyle m_{\bar{a},\bar{\pi}}^{2}$ $\displaystyle=$
$\displaystyle\frac{1}{2}\bigg{[}(m_{a}^{2}+m_{\pi}^{2})\mp$
$\displaystyle(m_{\pi}^{2}-m_{a}^{2})\sqrt{1+\frac{4\varepsilon^{2}\,\text{min}\\{m_{\pi}^{4},m_{a}^{4}\\}}{(m_{\pi}^{2}-m_{a}^{2})^{2}}}\bigg{]}~{},$

where

$\displaystyle m_{\pi}^{2}$ $\displaystyle=$ $\displaystyle
B_{0}(m_{u}+m_{d})~{},$ (21) $\displaystyle m_{a}^{2}$ $\displaystyle=$
$\displaystyle
m_{a_{0}}^{2}+4\,m_{\pi}^{2}\,\frac{m_{u}m_{d}}{(m_{u}+m_{d})^{2}}\left(\frac{f_{\pi}}{f}\right)^{2}c_{GG}^{2}~{}.$

In terms of these parameters, and working to linear order in $\varepsilon$ but
to all orders in $\alpha$ (to accommodate $m_{a}^{2}\simeq m_{\pi}^{2}$), we
find the physical axion and neutral pion couplings to nucleons $N=p,n$ and DM
$\chi$ to be

$\displaystyle\mathcal{L}$ $\displaystyle\supset$
$\displaystyle\bar{N}\,\gamma^{\mu}\gamma^{5}\left(\frac{g_{0}}{2}u^{(0)}_{\mu}+\frac{g_{A}}{2}\sigma^{3}u^{(1)}_{\mu}\right)N$
(22) $\displaystyle\ +\
\bar{\chi}\gamma^{\mu}\gamma^{5}u^{(\chi)}_{\mu}\chi~{},$

with the axial iso-scalar coupling $g_{A}=\Delta_{u}^{p}-\Delta_{d}^{p}\simeq
1.246$ Zyla  _et al._ (2020); Aoki  _et al._ (2021), the iso-vector coupling
$g_{0}=\Delta_{u}^{p}+\Delta_{d}^{p}\simeq 0.339$ Aoki  _et al._ (2021), and

$\displaystyle\begin{aligned}
u_{\mu}^{(\chi)}&=-\frac{1}{2f}\left[(c_{\alpha}+\widetilde{\varepsilon}\,s_{\alpha})\,\partial_{\mu}\bar{a}+(s_{\alpha}-\widetilde{\varepsilon}\,c_{\alpha})\,\partial_{\mu}\bar{\pi}\right]~{},\
\ \\\
u^{(0)}_{\mu}&=\frac{C_{0}}{2f}\left[(c_{\alpha}+\widetilde{\varepsilon}\,s_{\alpha})\,\partial_{\mu}\bar{a}+(s_{\alpha}-\widetilde{\varepsilon}\,c_{\alpha})\,\partial_{\mu}\bar{\pi}\right]~{},\\\
u^{(1)}_{\mu}&=\frac{1}{f_{\pi}}\left[-(s_{\alpha}-\widetilde{\varepsilon}\,c_{\alpha})\,\partial_{\mu}\bar{a}+(c_{\alpha}+\widetilde{\varepsilon}\,s_{\alpha})\,\partial_{\mu}\bar{\pi}\right]~{},\end{aligned}$
(23)

where $C_{0}=c_{uu}+c_{dd}+2c_{GG}$, and
$\widetilde{\varepsilon}=\varepsilon\;\Theta(m_{\pi}^{2}-m_{a}^{2})$.

The couplings in Eq. (23) imply that DM-nucleon interactions are mediated by a
combination of axion and neutral pion exchange. The leading Feynman diagrams
for spin-dependent $N+\chi\to N+\chi$ scattering are given in Fig. 3(a), and
the corresponding matrix element is

$\displaystyle\mathcal{M}$ $\displaystyle=$
$\displaystyle[\bar{u}_{N^{\prime}}i\gamma^{5}u_{N}]\,[\bar{u}_{\chi^{\prime}}i\gamma^{5}u_{\chi}]\,C_{a}^{N}\
,$ (24)

where

$\displaystyle
C_{a}^{N}=\frac{m_{\chi}m_{N}}{f^{2}}\,\frac{2\,c_{GG}}{q^{2}-m_{a}^{2}}\,{G}_{a}^{N}(q^{2})\
,$ (25)

with the form factor

$\displaystyle 4c_{GG}\,{G}_{a}^{N}(q^{2})$ (26) $\displaystyle=$
$\displaystyle\Bigg{(}g_{0}\,C_{0}\;\frac{q^{2}-c_{\alpha}^{2}m_{\pi}^{2}-s_{\alpha}^{2}m_{a}^{2}-2\widetilde{\varepsilon}s_{\alpha}c_{\alpha}\,(m_{\pi}^{2}-m_{a}^{2})}{q^{2}-m_{\pi}^{2}}$
$\displaystyle+\sigma^{3}_{N}\,g_{A}\left(\frac{2f}{f_{\pi}}\right)\,\left[s_{\alpha}c_{\alpha}-\widetilde{\varepsilon}(c_{\alpha}^{2}-s_{\alpha}^{2})\right]\,\frac{m_{\pi}^{2}-m_{a}^{2}}{q^{2}-m_{\pi}^{2}}\Bigg{)}$
$\displaystyle\ \to\ $ $\displaystyle
g_{0}C_{0}-\sigma^{3}_{N}\,g_{A}\left(\frac{2f}{f_{\pi}}\right)\varepsilon\,\frac{m_{\pi}^{2}}{q^{2}-m_{\pi}^{2}}\
,$

where $\sigma^{3}_{N}=\pm 1$ for $N=p,\,n$ and the last expression applies for
$|m_{a}^{2}-m_{\pi}^{2}|\gg\varepsilon\,m_{\pi}^{2}$. For $c_{uu},\,c_{dd}\to
0$, this form factor matches the leading CP-odd gluonic nucleon form factor
presented in Refs. Bishara  _et al._ (2017); Del Nobile (2021) and given in
Eq. (34) below.

In the non-relativistic limit, this matrix element leads to spin-dependent
scattering on nuclei with a per-nucleon scattering cross section of

$\displaystyle\frac{d\sigma_{\rm SD}}{dQ^{2}}$ $\displaystyle=$
$\displaystyle\frac{1}{16\pi}\,\frac{Q^{4}}{Q_{\rm
max}^{2}(m_{\chi}+m_{N})^{2}}\,|C_{a}^{N}|^{2}\ ,$ (27)

where $Q^{2}=|q^{2}|$ is the 3-momentum transfer, and $Q\leq
Q_{\text{max}}=2\,\mu_{N\chi}v_{\chi}$ for incoming non-relativistic DM
velocity $v_{\chi}$. A key feature of this result is that the total cross
section is suppressed by four powers of $v_{\chi}$ after integrating up to
$Q_{\rm max}^{2}$.

####  III.1.2 ALP Spin-Independent

The exchange of the radial mode $r$ leads to nuclear scattering independent of
spin and velocity, via the $t$-channel Feynman diagram in Fig. 3(b). Since the
radial mode is heavy, it can be safely integrated out based on the
interactions of Eq. (II.1.1) to yield an effective operator of the form
$(9\alpha_{s}/8\pi)\,C_{G}\,\bar{\chi}\chi\,GG$ with

$\displaystyle C_{G}=\frac{4\,c_{GG}}{27}\,\frac{m_{\chi}}{f^{2}m_{r}^{2}}\ .$
(28)

Matching this gluon operator to nucleons yields an effective SI-nucleon dark
matter cross section of

$\displaystyle\frac{d\sigma_{\rm SI}}{dQ^{2}}$ $\displaystyle=$
$\displaystyle\frac{1}{Q_{\rm
max}^{2}}\frac{\mu_{N\chi}^{2}}{\pi}\,|C^{N}_{r}|^{2}\ ,$ (29)

where $\mu_{N\chi}$ is the nucleon-DM reduced mass,
$C^{N}_{r}=C_{G}\,m_{N}f_{G}^{N}$, and $m_{N}f_{G}^{N}\simeq 0.85\,\text{GeV}$
is the CP-even gluon nucleon form factor Bishara  _et al._ (2017); Del Nobile
(2021). Integrating over $Q^{2}$ gives

$\displaystyle\sigma_{\rm SI}$ $\displaystyle\simeq$ $\displaystyle(2.0\times
10^{-46}\,\text{cm}^{2})\ \times\ $
$\displaystyle~{}~{}c_{GG}^{2}\left(\frac{\mu_{N\chi}\,m_{\chi}}{\text{GeV}^{2}}\right)^{2}\left(\frac{f}{m_{r}}\right)^{4}\left(\frac{100\,\text{GeV}}{f}\right)^{8}\
,$

Note that we have implicitly normalized the radial mode mass $m_{r}$ to $f$, a
benchmark we will use throughout this paper.

###  III.2 2HDM+$a$ Model Results

As in the ALP scenario, the 2HDM+$a$ scenario gives rise to both spin-
dependent and spin-independent nuclear scattering. The former arises from
exchange of the pseudoscalars and the latter via quantum loops, as depicted in
Fig. 3(a),(c). A detailed treatment follows.

####  III.2.1 2HDM+$a$ Spin-Dependent

The SD per-nucleon cross section with $t$-channel mediation by $a$ and $A$ as
in Fig. 3(a) may be written as

$\displaystyle\frac{d\sigma_{SD}}{dQ^{2}}=\frac{1}{16\pi}\frac{Q^{4}}{Q_{\rm
max}^{2}(m_{\chi}+m_{N})^{2}}\,|C_{\rm PS}^{N}|^{2}$ (30)

with

$\displaystyle C^{N}_{\rm
PS}=\sum_{q=u,d,s}C_{q}\,(G_{q}^{N}+G_{G}^{N})+\sum_{q=c,b,t}C_{q}\,G_{G}^{N}\
,$ (31)

where

$\displaystyle
C_{q}=\frac{m_{N}}{v}\bigg{(}\frac{\xi^{\chi}_{a}\xi_{a}^{q}}{m_{a}^{2}-q^{2}}+\frac{\xi^{\chi}_{A}\xi_{A}^{q}}{m_{A}^{2}-q^{2}}\bigg{)}\
,$ (32)

and

$\displaystyle G_{q}^{N}(q^{2})$ $\displaystyle=$
$\displaystyle\Delta_{q}^{N}-q^{2}\left(\frac{a_{q,\pi}^{N}}{q^{2}-m_{\pi}^{2}}+\frac{a_{q,\eta}^{N}}{q^{2}-m_{\eta}^{2}}\right)$
(33) $\displaystyle G_{G}^{N}(q^{2})$ $\displaystyle=$
$\displaystyle(-1)\\!\\!\sum_{q=u,d,s}\\!\frac{\bar{m}}{m_{q}}G_{q}^{N}(q^{2})\
.$ (34)

Here, $\bar{m}=(1/m_{u}+1/m_{d}+1/m_{s})^{-1}$, $m_{\pi}=135\,\text{MeV}$,
$m_{\eta}=548\,\text{GeV}$, and the remaining coefficients can be expressed in
terms of the quark spin fractions $\Delta_{q}^{N}$. We use the values from the
recent average of $2+1+1$ lattice results of Ref. Aoki  _et al._ (2021),

$\displaystyle\Delta_{u}^{p}=0.777\ ,\qquad\Delta_{d}^{p}=-0.483\
,\qquad\Delta_{s}^{p}=-0.053\ ,\qquad$ (35)

with $\Delta_{u}^{n}=\Delta_{d}^{p}$, $\Delta_{d}^{n}=\Delta_{u}^{p}$,
$\Delta_{s}^{n}=\Delta_{s}^{p}$, and $g_{A}=\Delta_{u}^{p}-\Delta_{d}^{p}$.
The pole coefficients are are Bishara  _et al._ (2017); Del Nobile (2021)

$\displaystyle a_{u,\pi}^{p}=-a_{d,\pi}^{p}=\frac{1}{2}g_{A}\
,~{}~{}a_{s,\pi}^{p}=0\ ,$

with $a_{u,\pi}^{n}=a_{d,\pi}^{p}$, $a_{d,\pi}^{n}=a_{u,\pi}^{p}$,
$a_{s,\pi}^{n}=a_{s,\pi}^{p}$, and

$\displaystyle
a_{u,\eta}^{p,n}=a_{d,\eta}^{p,n}=-\frac{1}{2}a_{s,\eta}^{p,n}=\frac{1}{6}(\Delta_{u}^{p}+\Delta_{d}^{p}-2\Delta_{s}^{p})\
.$

We note that these quantities have moderate uncertainties, but their effect is
significantly less than other uncertainties related to neutron star capture.

####  III.2.2 2HDM+$a$ Spin-Independent

In the minimal 2HDM+$a$ model presented above there is no tree-level
contribution to SI nucleon scattering. However, a velocity-independent SI
interaction is generated at loop level as illustrated in Fig. 3. These
contributions were studied in Refs. Ipek  _et al._ (2014); Arcadi  _et al._
(2018); Bell  _et al._ (2018b) and developed further in Refs. Abe  _et al._
(2019); Ertas and Kahlhoefer (2019).

The SI nucleon cross section can be written as

$\displaystyle\frac{d\sigma_{\text{SI}}}{dQ^{2}}=\frac{1}{Q_{\rm
max}^{2}}\frac{\mu^{2}_{N\chi}}{\pi}\left|C^{N}_{\rm loop}\right|^{2}$ (36)

where

$\displaystyle C^{N}_{\rm loop}\simeq
C_{G}\,m_{N}f^{N}_{G}+\sum_{q=u,d,s}C_{q}\,m_{N}f^{N}_{q}\ ,$

with the recent lattice determinations Alexandrou  _et al._ (2020)

$\displaystyle f_{u}^{p}=0.017\ ,\qquad f_{d}^{p}=0.026\ ,\qquad
f_{s}^{p}=0.049\ ,\qquad$ (37) $\displaystyle f_{u}^{n}=0.013\ ,\qquad
f_{d}^{n}=0.033\ ,\qquad f_{s}^{n}=0.049\ ,\qquad$

and $f_{G}^{N}=1-\sum_{q}f_{q}^{N}\simeq 0.91$.

Explicit expressions for the loop-induced $C_{q}$ and $C_{G}$ coefficients in
the 2HDM+$a$ are collected in Ref. Abe  _et al._ (2019). For the model
parameters considered in this work, the dominant loops for SI scattering
correspond to the lower two diagrams in Fig. 3. The $C_{q}$ coefficients in
Eq. (III.2.2) are generated mainly by triangle diagrams of the form shown in
Fig. 3 (c). The $C_{G}$ coefficient comes primarily from triangle and box
loops connecting to heavy ($c,b,t$) quarks. We consider $\tan\beta=3$, but at
larger values a contribution from two-loop diagrams involving gluons can come
to dominate due to a $\tan^{4}\\!\beta$ enhancement in loops involving bottom
quarks Abe  _et al._ (2019). Higher twist operators are also considered in
Ref. Abe  _et al._ (2019), but we find that their numerical contributions here
are very small.

Let us note that loop diagrams in the ALP benchmark do not contribute to the
SI cross section in this model to as much of an extent as in the 2HDM+$a$ .
There is no Higgs-ALP coupling generated at leading order Bauer  _et al._
(2021b), and this eliminates the triangle diagram. For the box diagrams, the
crucial difference is that the ALP couplings of Eqs. (II.1.1,7) involve
derivative operators rather than a true pseudoscalar-fermion interaction.
These are equivalent on-shell, but in the box loops the derivative couplings
add powers of loop momentum and make the result UV-dependent. The leading
contribution to SI scattering in the ALP benchmark is therefore from the
radial mode rather than loops.

Figure 4: Contours of constant NS luminosity, logarithmically equispaced, for
effective couplings $y_{\chi}^{2}y_{n}^{2}=10^{-14}=1$; here $L_{\rm max}$ is
the maximum NS luminosity that may be imparted by kinetic heating due to
incident dark matter flux. Equivalently, these are contours of
$y_{\chi}^{2}y_{n}^{2}$ corresponding to NS luminosity $L_{\rm max}$, with
values of $y_{\chi}^{2}y_{n}^{2}$ that increase as we go from inner regions of
the plot to the outer.

##  IV Neutron star heating

We turn next to study the heating of neutron stars by the infall and capture
of local $\chi$ dark matter based on the two mediator benchmark scenarios
introduced above. First we review general features of dark kinetic heating of
neutron stars, and then we specialize to the specific aspects of the capture
and heating when DM connects to the SM through a light pseudoscalar mediator.

###  IV.1 Review of Dark Kinetic Heating

For our treatment of DM capture in neutron stars we will consider the
following representative NS configuration:

$\displaystyle M_{\rm NS}=1.5\ M_{\odot},~{}~{}R_{\rm NS}=12.85\ {\rm km}\ ,$
(38)

giving an escape speed at the surface

$\displaystyle v_{\rm esc}=\sqrt{\frac{2GM_{\rm NS}}{R_{\rm NS}}}\simeq 0.59\
.$ (39)

This configuration is obtained from a Quark Meson Coupling (QMC) equation of
state (EOS) of matter near nuclear saturation densities in Ref. Anzuini  _et
al._ (2021), where it is also shown that near-identical results are obtained
for the Brussels-Montreal EOS “BSk24”. The self-interactions of the nuclear
medium in the NS core effectively modify the neutron mass,
$m_{n}\to\bar{m}_{n}$, and in our work we set

$\displaystyle\bar{m}_{n}=719.6\ {\rm MeV}\ ,$ (40)

obtained by volume-averaging the NS radius-dependent effective $m_{n}$ shown
in Ref. Anzuini  _et al._ (2021). We will use this representative NS
configuration to illustrate our main effects and to present the central values
of our parametric reaches while accounting for uncertainties in the (yet
unknown) mass and radius of candidate NSs.

To estimate the sensitivity of NS capture in the coupling-mass space of the
underlying DM theory, we first compute the fraction of the DM flux incident on
the NS that is gravitationally captured. For local DM density $\rho_{\chi}$
and average DM-NS relative speed $v_{\rm rel}$ (which at the solar position
are respectively taken as $0.4$ GeV/cm3 and 350 km/s Eilers  _et al._ (2019)),
the DM mass capture rate is given by Goldman and Nussinov (1989)

$\displaystyle\dot{M}=m_{\chi}C_{n\chi}=\rho_{\chi}v_{\rm rel}\times\pi b_{\rm
max}^{2}\times p_{v}\times p_{\sigma}\ ,$ (41)

where $b_{\rm max}=R_{\rm NS}(1+z)(v_{\rm esc}/v_{\rm rel})$ is the maximum
impact parameter of DM traversing the NS, with $1+z=(1-v_{\rm
esc}^{2})^{-1/2}$ a blueshift factor magnifying the NS radius to a distant
observer, and $p_{v}$ is the probability that a scattered DM particle loses
sufficient energy to be captured; in Appendix C we show that for pseudoscalar
mediation $p_{v}=1$ is an excellent approximation throughout the parameter
space we consider. The probability that incident DM is scattered is given by

$\displaystyle p_{\sigma}$ $\displaystyle=$ $\displaystyle
1-e^{-\tau}\simeq\tau$ (42) $\displaystyle=$
$\displaystyle\sigma_{n\chi}/\sigma_{\rm cap}\ ,$

where, for optical depth $\tau$, the approximate equality in the first line
holds in the optically thin limit considered here. The “capture cross section”
above which $\tau>1$ in the NS core is:

$\displaystyle\displaystyle\sigma_{\rm
cap}=\begin{cases}\sigma_{0}(\bar{m}_{n}/m_{\chi})\ &,\ \
m_{\chi}<\bar{m}_{n}\ ,\\\ \sigma_{0}\ &,\ \ \bar{m}_{n}\leq
m_{\chi}\leq\text{PeV}\ ,\\\ \sigma_{0}(m_{\chi}/{\rm PeV})\ &,\ \
m_{\chi}>{\rm PeV}\ ,\end{cases}$ (43)

where the NS geometric cross section $\sigma_{0}=\pi(\bar{m}_{n}/M_{\rm
NS})R_{\rm NS}^{2}\simeq 2.2\times 10^{-45}\,{\rm cm}^{2}$. The dependence on
$m_{\chi}$ in Eq. (43) is understood by considering the typical DM recoil
energy in the neutron rest frame:

$\displaystyle E_{\rm rec}$
$\displaystyle\simeq\frac{\bar{m}_{n}m_{\chi}^{2}\gamma^{2}v^{2}_{\rm
esc}}{(\bar{m}_{n}^{2}+m_{\chi}^{2}+2\gamma\bar{m}_{n}m_{\chi})}\ ,$ (44)

with $\gamma=1+z$. For $m_{\chi}\\!<\\!\bar{m}_{n}$, only a fraction $\simeq
3\Delta p/p_{F}$ of degenerate neutrons close enough to their Fermi surface
receive the typical momentum transfer $\Delta p=\sqrt{2\bar{m}_{n}E_{\rm
rec}}$ to scatter111As noted in Ref. Bramante  _et al._ (2022), for
$m_{\chi}\lesssim 35$ MeV nucleons in the NS might not scatter elastically due
to a superfluidity energy gap, leaving DM to capture via collective
excitations instead. A detailed treatment of this effect is beyond the scope
of this work, however, we note that in the DM mass range where it applies,
current limits on our mediator already outdo the future sensitivity of NS
heating: see Fig. 1. For a study of collective effects in NSs impacting light
DM capture in certain models, see Ref. DeRocco  _et al._ (2022). to above the
Fermi momentum $p_{F}\simeq 0.4~{}\text{GeV}$. This so-called Pauli blocking
gives $\sigma_{\rm cap}\propto E_{\rm rec}^{-1/2}\propto m_{\chi}^{-1}$. For
$\bar{m}_{n}\\!\leq\\!m_{\chi}\\!\leq\\!10^{6}~{}\\!\text{GeV}$, a single
scatter suffices for capture: $E_{\rm rec}\simeq\bar{m}_{n}v_{\rm
esc}^{2}\gamma^{2}>$ KEhalo, the DM halo kinetic energy. For
$m_{\chi}>\text{PeV}$, multiple scatters are required for capture, so that
approximately $\sigma_{\rm cap}\propto{\rm KE}_{\rm halo}/E_{\rm rec}\propto
m_{\chi}$.

Under equilibrium, the kinetic power of the infalling dark matter (the heating
rate) equals the rate at which photons are emitted from the NS surface (the
cooling rate, dominated by such photon emission for NSs older than $\sim$Myr
Page  _et al._ (2004); Yakovlev and Pethick (2004)). Assuming blackbody
radiation, the NS luminosity corresponding to a temperature $T$ (in the NS
frame) is then

$\displaystyle L=z\dot{M}=4\pi R_{\rm NS}^{2}T^{4}\ .$ (45)

This luminosity attains a maximum value $L_{\rm max}$ for unit capture
probabilities $p_{\sigma}$ and $p_{v}$.

Applying these general results to our representative NS configurationwith
parameters as in Eq. (38), we find the maximum luminosity

$\displaystyle L_{\rm max}=7.6\times 10^{24}~{}{\rm GeV/s}\ .$ (46)

This corresponds to a NS temperature seen by a distant observer
$\widetilde{T}=T/(1+z)$ of

$\displaystyle\widetilde{T}=1400~{}{\rm K}\ .$ (47)

Such a NS temperature is measurable within reasonable integration times at
current and imminent infrared telescope missions Baryakhtar  _et al._ (2017);
Raj  _et al._ (2018), in particular at the recently launched James Webb Space
Telescope (JWST)222See Ref. Chatterjee  _et al._ (2022) for an astronomical
study on observing NS heating at JWST. Gardner  _et al._ (2006), and the
forthcoming Extremely Large Telescope Neichel  _et al._ (2018) and Thirty
Meter Telescope Skidmore  _et al._ (2015).

Possible NS-reheating mechanisms that may compete with DM kinetic heating are:
(a) the accretion of interstellar material, likely to be deflected along the
NS’s magnetic field lines such that only a small polar region is heated,
distinguishable from all-surface thermal emission Treves  _et al._ (2000); (b)
rotochemical heating in the case of NSs with very small (sub-7 ms) initial
spin period for certain nucleon pairing models Hamaguchi  _et al._ (2019a).

While we focus on the heating of neutron stars from the transfer of kinetic
energy by captured dark matter in this work, which applies to any type of dark
matter candidate, we note that dark matter that annihilates within a NS could
raise its brightness further Kouvaris (2008); de Lavallaz and Fairbairn (2010)
and thereby improve telescope sensitivities Baryakhtar  _et al._ (2017); Raj
_et al._ (2018); Chatterjee  _et al._ (2022). This effect depends non-
trivially on whether or not the DM thermalizes with the NS over the stellar
lifetime Bertoni  _et al._ (2013); Garani  _et al._ (2021), and we defer its
study to a future analysis.

###  IV.2 Capture via (Pseudo)scalar Mediation

In Sec. III we computed per-nucleon scattering cross sections relevant for
direct detection experiments. To apply these to DM capture in neutron stars,
where the incoming DM and target nucleons are more energetic and can both be
semi-relativistic Joglekar  _et al._ (2020b); Bell  _et al._ (2021a); Anzuini
_et al._ (2021), we generalize our previous results in two ways. First, we
extend the nucleon matrix element coefficients $C_{i}^{N}$ of Eqs.
(26,28,32,III.2.2)by multiplying them by a dipole form factor Bishara  _et
al._ (2017); Del Nobile (2021),

$\displaystyle C_{i}^{N}(t)\ \to\
C_{i}^{N}(t)\times\frac{1}{(1-t/M_{D}^{2})^{2}}$ (48)

with $t=q^{2}=-Q^{2}$ and $M_{D}\simeq 1\,\text{GeV}$. This form factor
accounts for the inner structure of the nucleon probed at higher momenta and
corresponds to expanding beyond the leading non-trivial order in an expansion
in $t/m_{N}^{2}$ Bishara  _et al._ (2017); Del Nobile (2021). Our second
modification is to use fully relativistic expressions for the cross sections.
In the notation of Sec. III, the summed, squared, and averaged matrix elements
are

$\displaystyle|\mathcal{M}|^{2}=\left\\{\begin{array}[]{lcl}|C_{i}^{N}|^{2}\,(4m_{\chi}^{2}-t)\,(4m_{N}^{2}-t)&;&\text{scalar}\\\
&&\\\ |C_{i}^{N}|^{2}\,t^{2}&;&\text{pseudoscalar}\end{array}\right.$ (52)

where $t=q^{2}=-Q^{2}$ is the usual Mandelstam variable.

Consider now the scattering of DM with neutrons mediated by $t$-channel
pseudoscalar exchange. The dependence of the cross sections on momentum
transfer makes this effectively unobservable in direct detection, but we
demonstrate below that these channels can be dominant for DM capture in
neutron stars. To illustrate how pseudoscalar mediation can be enhanced in NS
capture with a focus on the impact of the mediator mass, we use a simplified
pseudoscalar model with $C_{\rm PS,a}^{N}\to y_{\chi}y_{n}/(m_{a}^{2}-t)$ for
effective DM and neutron couplings $y_{\chi}$ and $y_{n}$. In Fig. 4 we show
contours of constant NS luminosity, corresponding one-to-one to contours of
the capture probability $p_{\sigma}$, in the plane of $m_{\chi}$ versus
$m_{a}/m_{\chi}$ for the fixed coupling product
$y_{\chi}^{2}y_{n}^{2}=10^{-14}$. Here, the entire region below the contour
depicting $L_{\rm max}$ corresponds to $p_{\sigma}=1$. The scaling behaviour
of these contours can be understood as follows. For $m_{\chi}\ll m_{n}$, the
typical momentum transfer is $Q\simeq m_{\chi}v_{\rm esc}(1+z)\simeq
m_{\chi}$, so we have

$\displaystyle\sigma_{n\chi}|_{m_{\chi}\ll m_{n}}$ $\displaystyle\propto$
$\displaystyle\frac{m_{\chi}^{4}}{(m_{\chi}^{2}+m_{a}^{2})^{2}}~{}.$ (53)

Combining this with the fact that for $m_{\chi}<m_{n}$ Pauli-blocking sets
$p_{\sigma}\propto m_{\chi}^{-1}$ (Eqs. (42) and (43)), we obtain

$\displaystyle L|_{m_{\chi}\ll m_{n}}\propto
p_{\sigma}\propto\begin{cases}1/m_{\chi},~{}~{}~{}~{}m_{a}\ll m_{\chi},\\\
m_{\chi}^{3}/m_{a}^{4},~{}m_{a}\gg m_{\chi}\ .\end{cases}$ (54)

This is indeed the scaling we observe in the $m_{a}\ll m_{\chi}$ and $m_{a}\gg
m_{\chi}$ limits in Fig. 4.

For $m_{\chi}\gg m_{n}$, the typical momentum transfer is $Q\simeq m_{n}v_{\rm
esc}(1+z)\simeq m_{n}$, and hence

$\displaystyle\sigma_{n\chi}|_{m_{\chi}\gg m_{n}}$ $\displaystyle\propto$
$\displaystyle\frac{1}{m_{\chi}^{2}(m_{n}^{2}+m_{a}^{2})^{2}}\ .$ (55)

In this mass range there is no Pauli-blocking; hence we find these scalings in
Fig. 4:

$\displaystyle L|_{m_{\chi}\gg m_{n}}\propto
p_{\sigma}\propto\begin{cases}1/m_{\chi}^{2},~{}~{}~{}~{}m_{a}\ll m_{n},\\\
1/m_{\chi}^{2}m_{a}^{4},~{}m_{a}\gg m_{n}~{}.\end{cases}$ (56)

We will find the above features reflected in the NS heating sensitivities we
estimate below.

##  V Results

Having investigated our benchmark pseudoscalar-mediated DM models as well as
the general features of DM kinetic heating of neutron stars, we now present
our specific results in detail.

###  V.1 ALP Scenario

In Figure 1 we show in red and purple bands the sensitivity of NS heating to
the dimensionless parameter $c_{GG}\,m_{n}/f$ characterizing the ALP-nucleon
coupling strength as a function of the DM mass $m_{\chi}$. The solid curves
are the sensitivities to spin-dependent scattering on NS nucleons via ALP
exchange, whereas the dashed curves are to spin-independent scattering via
exchange of the radial mode $r$. The left panel corresponds to fixing the
mediator mass to $m_{a}=m_{\chi}/10$, representative of a “light mediator”
regime, whereas the right panel depicts $m_{a}=m_{\chi}$ to represent an
“intermediate mediator mass” regime.

The bands correspond to the reach over the model parameter space from the
measurement of a NS with luminosity $L_{\rm max}$ given in Eq. (46), where the
upper (lower) boundary is for a NS mass of $M_{\rm NS}=1.5\,M_{\odot}$
($2.16\,M_{\odot}$). In general, the precise mass of the NS under observation
is not known, with heavier NSs yielding a greater luminosity for a given
capture cross section. This is the largest astrophysical uncertainty in the
capture rate, spanning an order of magnitude and exceeding those coming from
other effects such as the unknown EoS of NS matter, the superfluid nature of
NS nucleons, etc. Berti  _et al._ (2022). We estimate that the capture cross
section varies by only a factor of 1.6 across NS mass-radius configurations
across the mass range $[1.5,\,2.16]\,M_{\odot}$ Anzuini  _et al._ (2021).
Therefore we may infer that, for an observed NS luminosity, there is an
uncertainty of $\mathcal{O}(10)$ in the minimum ALP-mediated DM cross section
deduced from it. In practice we fix this uncertainty to $\sigma_{\rm uncert}$
= 10, and as the ALP exchange ($r$ exchange) cross section $\propto f^{-4}$ (
$\propto f^{-8}$), the bands spread vertically over factors of $\sigma_{\rm
uncert}^{1/4}$ ($\sigma_{\rm uncert}^{1/8}$).

The solid curves are seen to change slope at two different DM masses: near
$m_{\chi}=$ 20 MeV in a sharp manner, and near $m_{\chi}=$ GeV more gently.
The first feature is due to cancellations between the iso-scalar and iso-
vector terms in the form factor in Eq. (26) for small momentum transfers. As
we raise $m_{\chi}$ above $m_{\pi}$ (while keeping $m_{\chi}<\bar{m}_{n}$),
Eq. (26) is dominated by the iso-scalar term, and we have $\sigma_{\rm
n\chi}\propto|\mathcal{M}|^{2}/\bar{m}_{n}^{2}\propto f^{-4}m_{\chi}^{2}$. In
this regime the NS capture cross section $\sigma_{\rm cap}\propto
m_{\chi}^{-1}$ (Eq. (43)), thus the constraint on $1/f\propto
m_{\chi}^{-3/4}$. As we transition to $m_{\chi}$ above $\bar{m}_{n}$, we have
$Q\simeq\bar{m}_{n}$, thus the spin-averaged $|\mathcal{M}|^{2}\propto
m_{\chi}^{2}/f^{4}/(\bar{m}_{n}^{2}+m_{a}^{2})^{2}$. For $m_{a}=m_{\chi}$, we
then have $\sigma_{\rm n\chi}\propto|\mathcal{M}|^{2}/m_{\chi}^{2}\propto
f^{-4}m_{\chi}^{-4}$. As the NS capture cross section here is
$m_{\chi}$-independent, the constraint on $1/f\sim m_{\chi}$. For
$m_{a}=m_{\chi}/10$ the turnover occurs at $m_{\chi}=$ a few
$\times\bar{m}_{n}$ as expected. We also see that ALP exchange dominates the
sensitivity to NS capture over exchange of the heavier $r$-mode, and that the
latter has a different slope due to the simple
$\mu_{n\chi}^{2}m_{\chi}^{2}/f^{8}$ dependence on the cross section.

Comparing across the two panels, we also see that near $m_{\chi}\simeq 10$
GeV, the constraints on $c_{GG}\,m_{n}/f$ are about 10 times weaker for
$m_{a}/m_{\chi}=1$ than for $m_{a}/m_{\chi}=1/10$. This is because in this
region the cross section goes as $f^{-4}m_{a}^{-4}$. However, as we go to
smaller $m_{\chi}$, and hence smaller $m_{a}$, the ALP propagator in Eq. (24)
is dominated by $Q^{2}$ and the cross section becomes insensitive to $m_{a}$,
resulting in similar limits. There is no change in the dashed curves across
the panels due to the $r$-exchange cross section being independent of $m_{a}$
(Eq. (29)). The brown region, which also does not change across the panels, is
the upper bound from direct detection searches at SuperCDMS Al-Bakry  _et al._
(2022), Xenon-1T Aprile  _et al._ (2019), DarkSide-50 Agnes  _et al._ (2018),
and PandaX-4T Meng  _et al._ (2021). These limits are obtained from the spin-
independent $r$-exchange; the spin-dependent ALP exchange is
$Q^{4}$-suppressed and yields much weaker limits. The direct detection limits
are seen to be generically weaker than the NS capture sensitivity due to their
large exclusion cross sections at sub-nuclear DM masses, where these searches
are limited by recoil thresholds.

The gray shaded regions in the backgrounds of these plots are excluded by
constraints on the ALP mediator from meson decay searches in beam dumps, and
from collider searches; see Appendix A for particulars. The beam dump searches
lose their sensitivity above $m_{a}$ of $\mathcal{O}$(100
MeV)$-\mathcal{O}$(GeV) due to kinematic limitation in producing mesons.
Likewise, for $m_{a}\gtrsim$ GeV a CMS search for the chromo-magnetic dipole
moment of the top quark Sirunyan  _et al._ (2019a) is limited by backgrounds,
yielding an upper bound of $c_{GG}\,m_{n}/f\lesssim 5\times 10^{-2}$. Thus NS
heating could emerge the sole probe in this mass regime for $c_{GG}\,m_{n}/f$
between $\sim 10^{-3}-0.05$, as seen in Fig. 1. A notable feature in these
limits is the islands of parameter space near $m_{a}=100$ MeV unconstrained by
beam dumps and colliders. We see that the upper island ranging between
$c_{GG}\,m_{n}/f\simeq 7\times 10^{-4}-3\times 10^{-2}$ will be entirely
covered by NS heating in the light mediator case, and mostly covered in the
equal-mass mediator case. The lower islands around $m_{a}\simeq 10^{-4}$ may
be partially covered by NS heating in the light mediator case; with improved
telescope sensitivities to lower NS luminosities, it is possible for some
coverage of these islands in both the light and equal-mass mediator scenarios.

###  V.2 2HDM+$a$ Scenario

In Figure 2 we show in red and purple bands the sensitivity of NS heating to
the $a$-$A$ mixing angle $\theta$ as a function of the DM mass $m_{\chi}$; we
cut off the plot at $\theta=\pi/4$ in order to keep the lighter eigenstate $a$
(heavier eigenstate $A$) mostly the interacting field $a_{0}$ ($A_{0}$) after
electroweak symmetry breaking. The left (right) panel corresponds to a “light”
(“intermediate”) mediator with $m_{a}/m_{\chi}=1/3$ ($m_{a}/m_{\chi}=1$). As
in the previous sub-section, the bands represent the uncertainty in NS mass-
radius configurations, taken to be $\sigma_{\rm uncert}=$ 10\. As the per-
nucleon scattering cross section $\sigma_{\rm n\chi}\propto\sin^{2}2\theta$
(Eq. (30)), the bands spread vertically over a factor of $\sqrt{\sigma_{\rm
uncert}}$ for small $\theta$. Further, as $\sigma_{\rm n\chi}\propto
m_{\chi}^{-6}$ for fixed $m_{a}/m_{\chi}$ ratios, and as the NS capture cross
section is $m_{\chi}$-independent in this regime (see Eq. (43)), the limits on
$\theta\propto m_{\chi}^{3}$ for small $\theta$, as seen in the figure.

The brown region in both panels is excluded by the most recent limits from
direct detection searches at PandaX-4T Meng  _et al._ (2021) and Xenon-1T
Aprile  _et al._ (2018), using the spin-independent loop-induced cross section
we derive in Sec. III.2.2. These correspond to exclusion cross sections that
range roughly between $4\times 10^{-47}$ cm2 – $4\times 10^{-46}\ {\rm
cm}^{2}$. The loops that result in SI scattering in direct detection result
also in SI scattering in NSs, but due to modest $Q$-dependence the NS heating
reach (for $\sigma_{\rm cap}\simeq 2\times 10^{-45}$ cm2) is weaker than the
direct detection limits, thus we do not present it here. We display contours
of constant (spin-independent) $\sigma_{\rm n\chi}=10^{-47},10^{-48},10^{-49}$
cm2, depicting the approximate march of progress expected of direct detection
experiments of the next two generations until the so-called neutrino floor is
reached Aalbers  _et al._ (2022). The gray regions in the background are
excluded by a host of collider and flavor data as detailed in Appendix B.

Comparing across the two panels, we find that the NS heating sensitivities on
$\theta$ are about ten times weaker in the right panel, where $a$ is thrice
heavier. This is just as expected from the scaling $\sigma_{\rm
n\chi}\propto\sin^{2}2\theta/m_{a}^{4}$ in Eq. (30), implying that the
exclusion on $\theta\propto m_{a}^{2}$ for small $\theta$. Our NS heating
reach is seen to be stronger than the $B\to K^{*}\mu\mu$ limits on the
mediator, applicable to $m_{a}\lesssim 4$ GeV, for $m_{\chi}/m_{a}=1/3$ in the
left panel. In this panel we also outdo collider limits on the mediator from
$\tau\tau\mu\mu$ searches for 12 GeV$\lesssim m_{\chi}\lesssim$ 25 GeV.
Whether we outdo ATLAS limits from undetected final states depends on the
unknown mass of the NS observed: a 1.5 $M_{\odot}$ NS does not, whereas
greater masses yielding greater NS luminosity (depicted by the red band) do.
In the right panel with $m_{\chi}/m_{a}=1$, our reach is stronger than these
limits by only a small margin near $m_{\chi}=$ 10 GeV for a 1.5 $M_{\odot}$
NS, while outdoing them for up to $m_{\chi}\simeq$ 16 GeV as we increase the
NS mass to 2.16 $M_{\odot}$.

As for SI direct detection, we observe that these limits are generally weaker
than those from flavor and colliders on the mediator, except for
$m_{\chi}(=m_{a})\gtrsim$ 62 GeV in the right panel, where there is a gap in
collider limits from a combination of kinematic limitations on Higgs decay and
lack of analysis on $bb\to A\to\tau\tau$ below $m_{\tau\tau}$ = 90 GeV. In any
case, due to weak dependence on DM mass in the loop-induced SI limits on
$\theta$ contrasted with the rapid climb ($\theta\propto m_{\chi}^{3}$) in the
NS heating reach, we are generally stronger than direct detection for small
$m_{\chi}$ in the mass range displayed, and weaker for higher $m_{\chi}$. It
is interesting to note that in low mass regions NS heating can reach smaller
couplings than future direct DM searches, even after the SI neutrino floor is
encountered in the latter. Thus in these regions ($m_{\chi}\lesssim$ 30 GeV
for $m_{\chi}/m_{a}=1/3$ and $m_{\chi}\lesssim$ 15 GeV for $m_{\chi}/m_{a}=1$)
NS heating may emerge a main, if not sole, means to test the 2HDM+$a$
framework.

##  VI Conclusions and Outlook

In this study we have shown that kinetic heating of neutron stars by dark
matter is a leading probe of dark matter candidates that connect to the SM
primarily through pseudoscalar mediators. We have done so within the context
of two well-motivated, self-consistent pseudoscalar mediator frameworks:
axion-like particles for sub-10 GeV DM masses, and the 2HDM+$a$ model for
larger DM masses. Our results are summarized in Figs. 1 and 2. These plots
demonstrate a strong complementarity between NS heating tests of these
scenarios and direct searches for DM and their associated mediators.

Our work motivates generalizations in a number of new directions:

  * •

It would be interesting to explore other pseudoscalar mediator frameworks
beyond the specific ALP and 2HDM$+a$ scenarios presented here. In particular,
while we have considered a KSVZ-type ALP that couples exclusively to gluons
(and DM) at the Peccei-Quinn breaking scale $f$, other ALP scenarios with
additional couplings to the SM would invoke further constraints Bauer  _et
al._ (2021c) and modify the effective nucleon coupling and hence the
sensitivity to NS heating. This could be realized in the context of DFSZ-type
ALPs Zhitnitsky (1980); Dine  _et al._ (1981) that could induce direct
couplings to quarks or leptons. In particular, primarily leptonic ALP
couplings could lead to DM capture via scattering with the large population of
leptons present in neutron stars Bell  _et al._ (2019); Garani  _et al._
(2019); Garani and Heeck (2019); Joglekar  _et al._ (2020a, b); Bell  _et al._
(2021b).

  * •

The motivation of this work is to identify and investigate self-consistent
scenarios unfriendly to direct detection but susceptible to NS capture.
Pseudoscalar mediation is one such important possibility. However, there are
other ones explored in the literature: models of inelastic DM (thermal
Higgsinos Baryakhtar  _et al._ (2017); Acevedo  _et al._ (2020) and other
electroweak multiplets Hamaguchi  _et al._ (2019b)), and muon-philic DM in a
gauged $L_{\mu}-L_{\tau}$ model Garani and Heeck (2019). Yet more
possibilities include DM in the keV-GeV mass range with small scattering cross
sections, DM with large nuclear cross sections shielded by the rock
overburden, composite DM with super-Planckian masses yielding too small a flux
for direct searches Bramante  _et al._ (2018, 2019a, 2019b); Adhikari  _et
al._ (2022); Carney  _et al._ (2022), models that lead to appreciable
clustering of DM so that the Earth encounters DM clumps too infrequently
Bramante  _et al._ (2022), and so on. We leave these avenues of exploration to
future authors.

  * •

Certain model-specific assumptions here may be relaxed, leading to a richer
phenomenology that we leave to future work. For instance, we chose the maximum
mediator mass to be the DM mass in order to avoid invisible decay modes,
however a recasting of collider and flavor constraints incorporating such
decays would be an interesting exercise. Similarly, whilst we had taken DM to
be asymmetric, turning on (even tiny) annihilation rates could not only
increase the resultant NS luminosity (see, e.g., Ref. Kouvaris (2008)), but
also invoke limits from indirect detection. An interesting computation in this
case would be of the timescale for DM thermalization with the NS. In
particular, in regions leading to NS capture via SI scattering via radial mode
mediation (regions above the dashed line in Fig. 1), repeated scatters with NS
nucleons would slow down DM particles and the $Q^{4}$-dependent ALP-mediated
SD cross section may get smaller than the SI cross section, leading to the
intriguing possibility that thermalization is driven initially by SD
scattering but later by SI scattering.

  * •

We have assumed here that either of our scenarios makes up the entire DM
population. Were it instead to constitute a fraction, the limits on the direct
detection cross section would weaken proportionally, but this does not
necessarily apply to NS heating. Rather, the DM-heated NS’ luminosity would be
proportionally smaller, and hence a detailed astronomical treatment such as in
Refs. Baryakhtar  _et al._ (2017); Raj  _et al._ (2018); Chatterjee  _et al._
(2022) must be carried out to determine the most optimal filter, integration
times, and so forth.

As with the overheating of neutron stars in complete models of exotic baryons
McKeen  _et al._ (2021a, b); Goldman  _et al._ (2022), their overheating by
the capture of dark matter in complete models of the mediator provides for
upcoming astronomical missions the added motivation of directly probing
fundamental physics.

## Acknowledgments

This work is supported in part by the Natural Sciences and Engineering
Research Council (NSERC) of Canada. TRIUMF receives federal funding via a
contribution agreement with the National Research Council (NRC) of Canada.

##  Appendix A Limits on the ALP Model

Figure 5: Limits on the ALP mediator in the plane of the effective ALP-nucleon
coupling versus the ALP mass, from beam dump, flavor, collider, and
astrophysical measurements. See text for further details.  Figure 6: Derived
limits on the 2HDM+$a$ model in terms of the light pseudoscalar mass $m_{a}$
and the pseudoscalar mixing angle $|\theta|$ from flavor tests, Higgs
measurements, and collider searches. The other 2HDM+$a$ model parameters are
kept fixed at the benchmark values discussed in the text with Type-II
couplings to fermions, $\tan\beta=3$, and
$m_{H}=m_{A}=m_{H^{+}}=600\,\text{GeV}$.

In this appendix we elaborate on the various limits on the ALP mediator as
shown in Fig. 5.

Collider and Beam Dump Limits

Many of the limits we display were described in Ref. Bauer  _et al._ (2021c)
in the space of $c_{GG}/f$ vs $m_{a}$; we reproduce them here in the
$c_{GG}\,m_{n}/f$ vs $m_{a}$ plane by making the substitution $c_{GG}\to
T_{2}(r)=1/2$. The key phenomenological divider here is the QCD scale $\sim
m_{\pi}$. A non-zero $c_{GG}$ induces sizeable ALP couplings to the photon via
the QCD anomaly, $c_{\gamma\gamma}\simeq-1.92c_{GG}$ Bauer  _et al._ (2021a),
so that for $m_{a}<m_{\pi}$ the diphoton mode $a\to\gamma\gamma$ dominates the
ALP decay width. For $m_{a}>m_{\pi}$ the branching to hadrons dominates, but
decays to muons $a\to\ell^{+}\ell^{-}$ at the sub-percent level turn out to be
important.

For $m_{a}<m_{\pi}$ the strongest terrestrial constraints come from searches
for $K^{+}\to\pi^{+}+X_{\rm inv}$ at NA62 Cortina Gil  _et al._ (2021), with
$X_{\rm inv}$ an invisible state. In this range of $m_{a}$ the
$a\to\gamma\gamma$ decay length exceeds the 10 m detector size and is thus
unseen Bauer  _et al._ (2021c). Limits from $K_{L}\to\pi^{0}X_{\rm inv}$ at
KOTO Ahn  _et al._ (2019) are weaker due to suppressed CP-conserving
couplings, and even weaker limits are set by $B\to K^{*}X_{\rm inv}$ at Belle
Grygier  _et al._ (2017) and $\pi^{+}\to ae^{+}\nu_{e}$ at PIENU Aguilar-
Arevalo  _et al._ (2021). For somewhat large $m_{a}$ and $c_{GG}$ the
$a\to\gamma\gamma$ mode becomes visible, and limits are obtained from
$K^{+}\to\pi^{+}\gamma\gamma$ and $K_{L}\to\pi^{0}\gamma\gamma$ at E949, NA48,
NA62, and KTeV Artamonov  _et al._ (2005); Lai  _et al._ (2002); Lazzeroni
_et al._ (2014); Abouzaid  _et al._ (2008).

For $2m_{\mu}<m_{a}<m_{B}$ limits from $B^{+}\to K^{+}a(\mu\mu)$ and $B^{0}\to
K^{*}a(\mu\mu)$ at LHCb Aaij  _et al._ (2015, 2017) apply, and for
$m_{a}\gtrsim$ GeV the strongest limits come from those placed by CMS Sirunyan
_et al._ (2019a) on the chromo-magnetic dipole moment of the top
$\hat{\mu}_{t}$, defined by the effective operator

$\displaystyle\mathcal{L}\in-\hat{\mu}_{t}\frac{g_{s}}{2m_{t}}\bar{t}\sigma^{\mu\nu}T^{a}tG^{a}_{\mu\nu}~{}.$
(57)

The ALP contributes to $\hat{\mu}_{t}$ at one-loop order via diagrams
proportional to $c^{2}_{tt}$ and $c_{tt}c_{GG}$, and via wavefunction
renormalization of the external gluon $\propto c^{2}_{GG}$ Bauer  _et al._
(2021c). For $\Lambda_{\rm QCD}\lesssim m_{a}<800$ MeV we have limits from the
visible mode $a\to\gamma\gamma$ at the CHARM beam dump Bergsma  _et al._
(1985), recast by F. Kling in Ref. Agrawal  _et al._ (2021). Near
$m_{a}=m_{\eta^{0}}\simeq 500$ MeV, due to strong ALP-$\eta^{0}$ mixing the
ALP decay length falls short of the detector distance, resulting in a gap in
the limits.

Supernova Limits

For small values of $c_{GG}\,m_{n}/f$ limits from data on supernovae come into
play. For $m_{a}$ up to around 200 MeV, the classic limit from the over-
cooling of SN 1897A by the application of the so-called Raffelt criterion is
shown in Fig. 5. Here we have used the limit derived in Ref. Ertas and
Kahlhoefer (2020) assuming ALPs are produced in the supernova via
bremsstrahlung in nucleon scattering due to non-zero $c_{GG}$. In principle,
ALPs are also produced by Primakoff conversion, $\gamma p\to ap$, however the
limits from this process assuming non-zero $c_{\gamma\gamma}$ (trivializing
all other couplings) turn out to be much weaker than those from nucleon
bremsstrahlung Ertas and Kahlhoefer (2020), therefore using just the latter is
a good approximation here. We also show limits from the non-observation at the
SMM satellite of gamma rays produced by long-lived ALPs decaying to photons
outside the proto-neutron star of SN 1987A, the so-called “ALP burst” Caputo
_et al._ (2022b). Both the cooling and ALP burst limits have a ceiling from
the fact that for sufficiently large couplings the ALP is trapped within the
proto-neutron star. Finally, we also display limits recently placed in Ref.
Caputo  _et al._ (2022a) by considering the boost in explosion energy from
$a\to\gamma\gamma$ in several low-energy supernovae observed. We note that
other astrophysical limits apply for coupling ranges below our description
above in the $m_{a}\lesssim 100$ MeV region Caputo  _et al._ (2022b), but we
do not display them here as the limits we show are already stronger than the
NS heating sensitivities in this mass range (see Fig. 1).

Bounds on the UV Model

Our ALP model also contains heavy vector-like mediator quarks $Q$ and the
radial mode $r(x)$. Both imply additional but model-dependent bounds on the
theory.

With the minimal quantum numbers $\mathbf{(3,1,0)}$ we have considered so far,
the $N_{Q}$ mediator quarks $Q$ are stable. This makes them dangerous if they
were fully thermalized in the early Universe. Such states can annihilate very
efficiently to produce small relic densities relative to dark matter for
$m_{Q}\lesssim 10\,\text{TeV}$ Kang  _et al._ (2008); Jacoby and Nussinov
(2007); De Luca  _et al._ (2018); Gross  _et al._ (2019). After QCD
confinement, the heavy quarks form bound color-neutral states with themselves
as well as with ordinary quarks. The analysis of Ref. De Luca  _et al._ (2018)
argues that rearrangement reactions convert the vast majority of $Q$-SM mixed
states to more deeply bound, neutral $QQQ$ and $Q\bar{Q}$ modes. Despite these
multiple depletion mechanisms, a small relic population of mixed $Q$-SM bound
states is left over. These carry non-integer electric charges and face
extremely strong constraints of $n_{Q\text{-SM}}/n_{B}\lesssim 10^{-21}$ from
searches for fractional charges in oil drop experiments using standard and
meteoritic materials Kim  _et al._ (2007); Perl  _et al._ (2009). Since we do
not address the cosmological history of dark matter or the mediators in this
work, we do not pursue these limits further and assume implicitly that no
significant relic density of $Q$s was created.

Direct collider searches for stable $Q$ quarks place bounds that are
independent of cosmology. Following $Q\bar{Q}$ pair production, the heavy
quarks would be expected to hadronize and lose energy as they pass through the
detector Fairbairn  _et al._ (2007), in analogy to long-lived $R$-hadrons
considered in the context of (very) split supersymmetry Wells (2003); Arkani-
Hamed and Dimopoulos (2005); Giudice and Romanino (2004). Generalizing the
calculations of Ref. Preuss and Valencia (2021) to color-triplet fermions, we
estimate that the LHC searches of Refs. Khachatryan  _et al._ (2016); Aaboud
_et al._ (2019) translate into $m_{Q}\gtrsim 1500\,\text{GeV}$. Note that
since $m_{Q}/(f/c_{GG})=y_{Q}\,N_{Q}/2\sqrt{2}$, this can be achieved for
$c_{GG}/f\lesssim 10^{-2}\,\text{GeV}^{-1}$.

Collider searches also test the ALP and its related radial mode. Both produce
jet final states but in different ways. For the ALP masses $m_{a}\leq
10\,\text{GeV}$ considered the dominant ALP decay mode is to dijets or
hadrons, and the signal becomes very challenging at hadron colliders due to
backgrounds. At $m_{a}=10\,\text{GeV}$, the leptophobic $Z^{\prime}$ search of
Ref. Sirunyan  _et al._ (2019b) translates into the bound $c_{GG}/f\lesssim
3\times 10^{-2}\,\text{GeV}^{-1}$, but we do not know of any useful bounds at
lower masses. The heavier radial mode with $m_{r}\sim f$ decays primarily via
$r\to aa$, leading to pairs of highly boosted dijets. Assuming that each of
these pairs is reconstructed as a single jet, the dijet search of Ref.
Sirunyan  _et al._ (2018a) gives no bound for $N_{Q}=1$ and $c_{GG}/f\lesssim
2\times 10^{-2}\,\text{GeV}^{-1}$ for $N_{Q}=10$ and $m_{r}=f$. We estimate
that monojet-like searches for $r\to\chi\bar{\chi}$ decays such as Ref. Aad
_et al._ (2021a) would give a similar sensitivity.

Motivated by the very strong bounds on stable mediator quarks, we also
consider a second scenario in which they also carry hypercharge $Y=1/3$. This
permits a small mass mixing with $d_{R}$-type quarks that allows their decays.
Hypercharging the mediators also leads to couplings of the ALP and radial
scalars to photons of the form

$\displaystyle\mathcal{L}_{eff}$ $\displaystyle\supset$
$\displaystyle\frac{\alpha}{4\pi}\,\bigg{(}2N_{c}Y^{2}\bigg{)}\,c_{GG}\frac{a}{f}\,F_{\mu\nu}\widetilde{F}^{\mu\nu}$
$\displaystyle+\frac{\alpha}{2\pi}\,\bigg{(}\frac{4N_{c}Y^{2}}{3}\bigg{)}\,c_{GG}\,\frac{r}{f}\,F_{\mu\nu}{F}^{\mu\nu}\
.$

For the ALP, this corresponds to $g_{a\gamma\gamma}=2\alpha\,c_{GG}/(3\pi
f)\simeq 1.7\times 10^{-3}\,c_{GG}/f$ in the standard notation Di Luzio  _et
al._ (2020).

On their own, these couplings to photons (or hypercharge vector bosons more
generally) do not give stronger constraints on the ALP than the gluon
couplings we have considered Bauer  _et al._ (2021b). When both photon and
gluon interactions are present, as we have here, bounds that rely on on-shell
$a\to\gamma\gamma$ are weakened further by competing decays to hadrons while
collider limits can become much stronger Mariotti  _et al._ (2018); Cid Vidal
_et al._ (2019); Knapen  _et al._ (2022). Extrapolating the recent LHC boosted
diphoton search of Ref. ATL (2022) to this scenario, we find $c_{GG}/f\lesssim
5\times 10^{-4}\,\text{GeV}^{-1}$ for $m_{a}=10\,\text{GeV}$. Unfortunately,
this study does not extend to lower ALP masses and therefore we do not obtain
a related limit in the $m_{a}\in[2,10)\,\text{GeV}$ mass range of primary
interest. The data scouting method proposed in Ref. Knapen  _et al._ (2022) is
projected to have sensitivity to $c_{GG}/f\sim 2\times
10^{-3}\,\text{GeV}^{-1}$. We have also studied the impact of the induced
couplings of the radial mode to photons of Eq. (A), but we do not find any
limits that are stronger than the dijet bounds discussed already because of
the dominance of $r\to aa$ decay mode together with the very small branching
ratio for $a\to\gamma\gamma$. Finally, we note that direct searches for pair-
produced heavy quarks give bounds that are similar to or slightly weaker than
$m_{Q}\gtrsim 1500\,\text{GeV}$ Aaboud  _et al._ (2018); ATL (2021a); CMS
(2022).

##  Appendix B Limits on the 2HDM+$a$ Model

The 2HDM+$a$ model consists of a 2HDM sector together with an additional
singlet pseudoscalar $a_{0}$. This latter state mixes with the $A_{0}$
pseudoscalar from the 2HDM to produce two physical CP-odd mass eigenstates $a$
and $A$ with mixing angle $\theta$. In our dark matter analysis we concentrate
primarily on the limit of $m_{a}\ll m_{A}$ with $|\theta|\ll 1$ such that the
lighter state is mostly singlet. As a benchmark for our study, we fix
$m_{A}=m_{H}=m_{H^{+}}=600\,\text{GeV}$, $\tan\beta=3$, $\lambda_{3}\simeq
m_{h}^{2}/v$, and $\cos(\beta-\alpha)=0$ corresponding to the so-called
alignment/decoupling limit Gunion  _et al._ (2000); Djouadi (2008). The dark
matter dynamics then depend on $m_{a}$ and $\theta$, which we scan over.

These benchmark parameters fix the properties of most of the 2HDM sector, and
they allow this sector to be consistent with direct and indirect bounds on it.
For $\tan\beta=3$, the heavy Higgs masses are large enough to be consistent
with flavour bounds from $B\to X_{s}\gamma$ Misiak and Steinhauser (2017) and
$B_{s}\to\mu\mu$ Abe  _et al._ (2019, 2020), as well as collider searches such
as $H^{+}\to t\bar{b}$ Aad  _et al._ (2021b) and $H/A\to\tau\bar{\tau}$
Sirunyan  _et al._ (2018b); Aad  _et al._ (2020a). Furthermore, taking
$m_{H}=m_{H^{+}}$ implies a residual custodial symmetry that allows the model
to reproduce the precision electroweak predictions of the SM to a very good
approximation Bauer  _et al._ (2017b). Similarly, the alignment limit yields
an SM-like Higgs boson with production and decay rates effectively identical
to the SM Gunion  _et al._ (2000); Djouadi (2008).

While the 2HDM sector is largely safe in this benchmark, the presence of a
lighter pseudoscalar implies bounds on the mixing angle $|\theta|$ for a given
mass $m_{a}$. As mentioned in Sec. II.2.2, flavor bounds can be important in
this benchmark as well for $m_{a}\lesssim 10\,\text{GeV}$. In our region of
parameter space, these are dominated by the flavor-changing coupling of $a$ to
a bottom and strange quark induced at one-loop Arcadi  _et al._ (2018). Using
the results of Arnan  _et al._ (2017), we compute the rate for
$B^{\pm}\to{K^{\pm}}^{(\ast)}(a\to\mu^{+}\mu^{-})$ and compare to the LHCb
search for promptly-decaying bosons in this channel Aaij  _et al._ (2015). In
addition, there is a subdominant constraint that comes from requiring that the
branching fraction for $B_{s}\to\mu^{+}\mu^{-}$ agree with the world average
of $(3.01\pm 0.35)\times 10^{-9}$ Zyla  _et al._ (2020). The resulting limits
on $\theta$ from these processes are shown in Fig. 6.

Pseudoscalar masses $m_{a}\lesssim m_{h}/2$ are also strongly constrained by
searches for rare SM Higgs decays. The decay width for $h\to aa$ in our
benchmark is Ipek  _et al._ (2014); Bauer  _et al._ (2017b); Abe  _et al._
(2019)

$\displaystyle\Gamma(h\to
aa)=\frac{\lambda_{haa}^{2}}{32\pi}\,m_{h}\,\sqrt{1-(2m_{a}/m_{h})^{2}}\ ,$
(59)

where

$\displaystyle\lambda_{haa}=\frac{(m_{A}^{2}-m_{a}^{2})}{2m_{h}\,v}\,\sin^{2}2\theta\
.$ (60)

With $m_{A}=600\,\text{GeV}$ this width exceeds the total SM Higgs width for
$|\theta|\gtrsim 0.05$ and thus Higgs measurements imply bounds on the mixing
angle near this level. Since we focus on $m_{a}\leq m_{\chi}$, the $a$
pseudoscalar decays visibly through its mixing with the 2HDM pseudoscalar $A$.
Expressions for the widths of the $a$ state are collected in Refs. Bauer  _et
al._ (2017b); Haisch  _et al._ (2018). Exotic decays of the SM Higgs have been
searched for the by LHC collaborations Cepeda  _et al._ (2021), and in Fig. 6
we show estimated limits from ATLAS searches for $h\to aa\to bb\mu\mu$ Aad
_et al._ (2022) and $h\to aa\to\tau\tau\mu\mu$ Aad  _et al._ (2015) as well as
CMS searches for $h\to aa\to bb\tau\tau$ Sirunyan  _et al._ (2018c) and $h\to
aa\to 4\mu$ Sirunyan  _et al._ (2019c). Non-standard Higgs decays are also
constrained by Higgs rate measurements if they do not contribute to the
standard SM Higgs search channels. Global fits to ATLAS Higgs measurements put
a limit on _undetected_ Higgs decays of BR${}_{\text{und}}<0.15$ Aad  _et al._
(2020b); ATL (2021b). We also show the corresponding limit on
$|\theta|\lesssim 0.04$ in Fig. 6 under the assumptions of unmodified SM Higgs
production and that $h\to aa$ decays do not significantly pollute SM Higgs
rates.

Direct searches at the LHC become more important for pseudoscalar masses above
the Higgs decay threshold. Collider bounds on the 2HDM+$a$ model have been
studied extensively for $m_{a}>2\,m_{\chi}$ with dominant $a\to\chi\bar{\chi}$
decays No (2016); Arcadi  _et al._ (2018); Haisch  _et al._ (2018); Abe  _et
al._ (2020), but they have not been investigated in as much detail when this
channel is closed and the pseudoscalar decays visibly. We find that LHC
searches for $H/A\to\tau\tau$ can be reinterpreted to give the strongest
limits on this scenario. For our benchmark of $\tan\beta=3$, the most
important production channel is in association with bottom quarks. We estimate
the production cross section by rescaling the $bb\phi$ cross sections
tabulated in Ref. de Florian  _et al._ (2016) by the appropriate factor of
$(\tan\beta\,\sin\theta)^{2}$, and derive bounds on $|\theta|$ by comparing to
the limits obtained from $\phi\to\tau\tau$ searches at ATLAS Aad  _et al._
(2020a) and CMS Sirunyan  _et al._ (2018b). The exclusions we find based on
these studies are shown in Fig. 6.

Figure 7: Top: Event distributions with respect to DM-nucleon scattering angle
variable $z\equiv\cos^{2}\theta_{\rm lab}$, bottom: probability of not
capturing in the NS following DM-nucleon scattering, left: pseudoscalar
mediation, right: real scalar mediation. Due to the $q^{4}$ dependence in the
cross section, the case of pseudoscalar mediators does not prefer forward
scatters with low energy transfers, unlike light scalar mediators. This
results in almost all scatters resulting in capture. See Appendix C for
further details.

##  Appendix C Light Mediator Kinematics
and NS Capture

In this appendix we show that the fraction of DM scattering events in the NS
that do not result in capture is highly negligible in our scenario, i.e.
$p_{v}\to 1$ in Eq. (41). The case of $p_{v}\ll 1$ could potentially occur if
the mediator is much lighter than the momentum transfer, favoring forward
scatters with small energy losses insufficient for capture. This can be seen
in the normalized event distribution with respect to the scattering angle
variable $z\equiv\cos^{2}\theta_{\rm lab}=(1-\cos\theta_{\rm CM})/2$,

$\displaystyle s(z)\equiv\frac{1}{\sigma}\frac{d\sigma}{dz}~{},$ (61)

plotted in the top panels of Fig. 7 with the left (right) panel corresponding
to mediation by a CP-odd (CP-even) scalar of mass $m_{\phi}$. Here we show
distributions for various values of $m_{\phi}/(2\mu_{n\chi})$ scanning between
the “light” and “heavy” mediator regimes, since the momentum transfer $Q\simeq
2\mu_{n\chi}\cos\theta_{\rm lab}$. We find agreement for the CP-even scalar
case with Ref. Dasgupta  _et al._ (2020).

We see that real scalar mediation is biased toward forward events with small
$z$ as the mediator gets lighter, as expected for
$d\sigma/dz\propto(Q^{2}+m_{\phi}^{2})^{-2}$. However, for pseudoscalar
mediation we see that events are biased toward larger $z$ as the mediator gets
lighter, with $s(z)$ getting flatter. This is because here $d\sigma/dz\propto
Q^{4}/(Q^{2}+m_{\phi}^{2})^{2}$; as $m_{\phi}$ is decreased, backward events
are more favored, and the sensitivity to $z$ is lost for $m_{\phi}\ll q$.

Capture in the NS is guaranteed if scattering on its constituents depletes DM
of its halo kinetic energy. Thus the probability of post-scattering capture is
defined by Dasgupta  _et al._ (2020)

$\displaystyle p_{v}$ $\displaystyle=$ $\displaystyle\int_{z_{\rm min}}^{1}dz\
s(z)~{},$ $\displaystyle z_{\rm min}$ $\displaystyle=$
$\displaystyle\frac{(m_{\chi}+m_{n})^{2}}{4m_{\chi}m_{n}}\frac{v_{\rm
halo}^{2}}{v_{\rm halo}^{2}+v_{\rm esc}^{2}}~{},$ (62)

where we set $v_{\rm halo}=10^{-3}c$ and $v_{\rm esc}=0.59c$ in practice. In
the bottom panels of Fig. 7 we plot the probability of not capturing,
$1-p_{v}$, for both the CP-odd and CP-even mediators in the light regime. As a
result of events weighted toward $z=1$, we find $1-p_{v}$ negligible, i.e.
$p_{v}=1$ is an excellent approximation in our treatment. This is not the case
for CP-even scalar mediators, where we see that it is possible for an
$\mathcal{O}(1)$ fraction of scattering events to evade capture for
$m_{\phi}/(2\mu_{\rm n\chi})\leq 0.01$, as also shown in Ref. Dasgupta  _et
al._ (2020).

## References

  * Goodman and Witten (1985) M. W. Goodman and E. Witten, Phys. Rev. D 31, 3059 (1985).
  * Jungman  _et al._ (1996) G. Jungman, M. Kamionkowski, and K. Griest, Phys. Rept. 267, 195 (1996), arXiv:hep-ph/9506380 .
  * Akerib  _et al._ (2017) D. S. Akerib _et al._ (LUX), Phys. Rev. Lett. 118, 021303 (2017), arXiv:1608.07648 [astro-ph.CO] .
  * Aprile  _et al._ (2018) E. Aprile _et al._ (XENON), Phys. Rev. Lett. 121, 111302 (2018), arXiv:1805.12562 [astro-ph.CO] .
  * Meng  _et al._ (2021) Y. Meng _et al._ (PandaX-4T), Phys. Rev. Lett. 127, 261802 (2021), arXiv:2107.13438 [hep-ex] .
  * Mitridate  _et al._ (2022) A. Mitridate, T. Trickle, Z. Zhang, and K. M. Zurek, in  _2022 Snowmass Summer Study_ (2022) arXiv:2203.07492 [hep-ph] .
  * Adhikari  _et al._ (2022) P. Adhikari _et al._ ((DEAP Collaboration)‡, DEAP), Phys. Rev. Lett. 128, 011801 (2022), arXiv:2108.09405 [astro-ph.CO] .
  * Carney  _et al._ (2022) D. Carney _et al._ , (2022), arXiv:2203.06508 [hep-ph] .
  * Kahn and Lin (2022) Y. Kahn and T. Lin, Rept. Prog. Phys. 85, 066901 (2022), arXiv:2108.03239 [hep-ph] .
  * Tucker-Smith and Weiner (2001) D. Tucker-Smith and N. Weiner, Phys. Rev. D 64, 043502 (2001), arXiv:hep-ph/0101138 .
  * Tucker-Smith and Weiner (2005) D. Tucker-Smith and N. Weiner, Phys. Rev. D 72, 063509 (2005), arXiv:hep-ph/0402065 .
  * Fan  _et al._ (2010) J. Fan, M. Reece, and L.-T. Wang, JCAP 11, 042 (2010), arXiv:1008.1591 [hep-ph] .
  * Freytsis and Ligeti (2011) M. Freytsis and Z. Ligeti, Phys. Rev. D83, 115009 (2011), arXiv:1012.5317 [hep-ph] .
  * Fitzpatrick  _et al._ (2013) A. L. Fitzpatrick, W. Haxton, E. Katz, N. Lubbers, and Y. Xu, JCAP 1302, 004 (2013), arXiv:1203.3542 [hep-ph] .
  * Dienes  _et al._ (2014) K. R. Dienes, J. Kumar, B. Thomas, and D. Yaylali, Phys. Rev. D90, 015012 (2014), arXiv:1312.7772 [hep-ph] .
  * Tinyakov  _et al._ (2021) P. Tinyakov, M. Pshirkov, and S. Popov, Universe 7, 401 (2021), arXiv:2110.12298 [astro-ph.HE] .
  * Berti  _et al._ (2022) E. Berti _et al._ , in  _2022 Snowmass Summer Study_ (2022) arXiv:2203.07984 [hep-ph] .
  * Goldman and Nussinov (1989) I. Goldman and S. Nussinov, Phys. Rev. D40, 3221 (1989).
  * Baryakhtar  _et al._ (2017) M. Baryakhtar, J. Bramante, S. W. Li, T. Linden, and N. Raj, Phys. Rev. Lett. 119, 131801 (2017), arXiv:1704.01577 [hep-ph] .
  * Bramante  _et al._ (2017) J. Bramante, A. Delgado, and A. Martin, Phys. Rev. D 96, 063002 (2017), arXiv:1703.04043 [hep-ph] .
  * Raj  _et al._ (2018) N. Raj, P. Tanedo, and H.-B. Yu, Phys. Rev. D 97, 043006 (2018), arXiv:1707.09442 [hep-ph] .
  * Acevedo  _et al._ (2020) J. F. Acevedo, J. Bramante, R. K. Leane, and N. Raj, JCAP 03, 038 (2020), arXiv:1911.06334 [hep-ph] .
  * Chen and Lin (2018) C.-S. Chen and Y.-H. Lin, JHEP 08, 069 (2018), arXiv:1804.03409 [hep-ph] .
  * Bell  _et al._ (2018a) N. F. Bell, G. Busoni, and S. Robles, JCAP 1809, 018 (2018a), arXiv:1807.02840 [hep-ph] .
  * Garani  _et al._ (2019) R. Garani, Y. Genolini, and T. Hambye, JCAP 1905, 035 (2019), arXiv:1812.08773 [hep-ph] .
  * Camargo  _et al._ (2019) D. A. Camargo, F. S. Queiroz, and R. Sturani, JCAP 1909, 051 (2019), arXiv:1901.05474 [hep-ph] .
  * Hamaguchi  _et al._ (2019a) K. Hamaguchi, N. Nagata, and K. Yanagi, Phys. Lett. B795, 484 (2019a), arXiv:1905.02991 [hep-ph] .
  * Keung  _et al._ (2020) W.-Y. Keung, D. Marfatia, and P.-Y. Tseng, JHEP 07, 181 (2020), arXiv:2001.09140 [hep-ph] .
  * Bell  _et al._ (2020) N. F. Bell, G. Busoni, S. Robles, and M. Virgato, JCAP 09, 028 (2020), arXiv:2004.14888 [hep-ph] .
  * Dasgupta  _et al._ (2020) B. Dasgupta, A. Gupta, and A. Ray, JCAP 10, 023 (2020), arXiv:2006.10773 [hep-ph] .
  * Garani  _et al._ (2021) R. Garani, A. Gupta, and N. Raj, Phys. Rev. D 103, 043019 (2021), arXiv:2009.10728 [hep-ph] .
  * Maity and Queiroz (2021) T. N. Maity and F. S. Queiroz, Phys. Rev. D 104, 083019 (2021), arXiv:2104.02700 [hep-ph] .
  * Bell  _et al._ (2021a) N. F. Bell, G. Busoni, T. F. Motta, S. Robles, A. W. Thomas, and M. Virgato, Phys. Rev. Lett. 127, 111803 (2021a), arXiv:2012.08918 [hep-ph] .
  * Zeng  _et al._ (2022) Y.-P. Zeng, X. Xiao, and W. Wang, Phys. Lett. B 824, 136822 (2022), arXiv:2108.11381 [hep-ph] .
  * Anzuini  _et al._ (2021) F. Anzuini, N. F. Bell, G. Busoni, T. F. Motta, S. Robles, A. W. Thomas, and M. Virgato, JCAP 11, 056 (2021), arXiv:2108.02525 [hep-ph] .
  * Bell  _et al._ (2019) N. F. Bell, G. Busoni, and S. Robles, JCAP 1906, 054 (2019), arXiv:1904.09803 [hep-ph] .
  * Garani and Heeck (2019) R. Garani and J. Heeck, Phys. Rev. D100, 035039 (2019), arXiv:1906.10145 [hep-ph] .
  * Joglekar  _et al._ (2020a) A. Joglekar, N. Raj, P. Tanedo, and H.-B. Yu, Phys. Lett. B, 135767 (2020a), arXiv:1911.13293 [hep-ph] .
  * Joglekar  _et al._ (2020b) A. Joglekar, N. Raj, P. Tanedo, and H.-B. Yu, Phys. Rev. D 102, 123002 (2020b), arXiv:2004.09539 [hep-ph] .
  * Bell  _et al._ (2021b) N. F. Bell, G. Busoni, S. Robles, and M. Virgato, JCAP 03, 086 (2021b), arXiv:2010.13257 [hep-ph] .
  * Fujiwara  _et al._ (2022) M. Fujiwara, K. Hamaguchi, N. Nagata, and J. Zheng, (2022), arXiv:2204.02238 [hep-ph] .
  * Hamaguchi  _et al._ (2022) K. Hamaguchi, N. Nagata, and M. E. Ramirez-Quezada, (2022), arXiv:2204.02413 [hep-ph] .
  * Nussinov (1985) S. Nussinov, Phys. Lett. B 165, 55 (1985).
  * Barr  _et al._ (1990) S. M. Barr, R. S. Chivukula, and E. Farhi, Phys. Lett. B 241, 387 (1990).
  * Kaplan  _et al._ (2009) D. E. Kaplan, M. A. Luty, and K. M. Zurek, Phys. Rev. D 79, 115016 (2009), arXiv:0901.4117 [hep-ph] .
  * Agrawal  _et al._ (2020) P. Agrawal, A. Parikh, and M. Reece, JHEP 10, 191 (2020), arXiv:2003.00021 [hep-ph] .
  * Peccei and Quinn (1977a) R. D. Peccei and H. R. Quinn, Phys. Rev. Lett. 38, 1440 (1977a).
  * Peccei and Quinn (1977b) R. D. Peccei and H. R. Quinn, Phys. Rev. D 16, 1791 (1977b).
  * Weinberg (1978) S. Weinberg, Phys. Rev. Lett. 40, 223 (1978).
  * Wilczek (1978) F. Wilczek, Phys. Rev. Lett. 40, 279 (1978).
  * Di Luzio  _et al._ (2020) L. Di Luzio, M. Giannotti, E. Nardi, and L. Visinelli, Phys. Rept. 870, 1 (2020), arXiv:2003.01100 [hep-ph] .
  * Kim (1979) J. E. Kim, Phys. Rev. Lett. 43, 103 (1979).
  * Shifman  _et al._ (1980) M. A. Shifman, A. I. Vainshtein, and V. I. Zakharov, Nucl. Phys. B 166, 493 (1980).
  * Bauer  _et al._ (2017a) M. Bauer, M. Neubert, and A. Thamm, JHEP 12, 044 (2017a), arXiv:1708.00443 [hep-ph] .
  * Bauer  _et al._ (2021a) M. Bauer, M. Neubert, S. Renner, M. Schnubel, and A. Thamm, JHEP 04, 063 (2021a), arXiv:2012.12272 [hep-ph] .
  * Bauer  _et al._ (2021b) M. Bauer, M. Neubert, S. Renner, M. Schnubel, and A. Thamm, Phys. Rev. Lett. 127, 081803 (2021b), arXiv:2102.13112 [hep-ph] .
  * Bauer  _et al._ (2021c) M. Bauer, M. Neubert, S. Renner, M. Schnubel, and A. Thamm, (2021c), arXiv:2110.10698 [hep-ph] .
  * Ertas and Kahlhoefer (2020) F. Ertas and F. Kahlhoefer, JHEP 07, 050 (2020), arXiv:2004.01193 [hep-ph] .
  * Caputo  _et al._ (2022a) A. Caputo, H.-T. Janka, G. Raffelt, and E. Vitagliano, Phys. Rev. Lett. 128, 221103 (2022a), arXiv:2201.09890 [astro-ph.HE] .
  * Caputo  _et al._ (2022b) A. Caputo, G. Raffelt, and E. Vitagliano, Phys. Rev. D 105, 035022 (2022b), arXiv:2109.03244 [hep-ph] .
  * Ipek  _et al._ (2014) S. Ipek, D. McKeen, and A. E. Nelson, Phys. Rev. D 90, 055021 (2014), arXiv:1404.3716 [hep-ph] .
  * Gunion  _et al._ (2000) J. F. Gunion, H. E. Haber, G. L. Kane, and S. Dawson,  _The Higgs Hunter’s Guide_ , Vol. 80 (2000).
  * Djouadi (2008) A. Djouadi, Phys. Rept. 459, 1 (2008), arXiv:hep-ph/0503173 .
  * Kling  _et al._ (2016) F. Kling, J. M. No, and S. Su, JHEP 09, 093 (2016), arXiv:1604.01406 [hep-ph] .
  * Arcadi  _et al._ (2018) G. Arcadi, M. Lindner, F. S. Queiroz, W. Rodejohann, and S. Vogl, JCAP 03, 042 (2018), arXiv:1711.02110 [hep-ph] .
  * Arnan  _et al._ (2017) P. Arnan, D. Bečirević, F. Mescia, and O. Sumensari, Eur. Phys. J. C 77, 796 (2017), arXiv:1703.03426 [hep-ph] .
  * Aaij  _et al._ (2015) R. Aaij _et al._ (LHCb), Phys. Rev. Lett. 115, 161802 (2015), arXiv:1508.04094 [hep-ex] .
  * Blinov  _et al._ (2022) N. Blinov, E. Kowalczyk, and M. Wynne, JHEP 02, 036 (2022), arXiv:2112.09814 [hep-ph] .
  * Zyla  _et al._ (2020) P. A. Zyla _et al._ (Particle Data Group), PTEP 2020, 083C01 (2020).
  * Aoki  _et al._ (2021) Y. Aoki _et al._ , (2021), arXiv:2111.09849 [hep-lat] .
  * Bishara  _et al._ (2017) F. Bishara, J. Brod, B. Grinstein, and J. Zupan, JHEP 11, 059 (2017), arXiv:1707.06998 [hep-ph] .
  * Del Nobile (2021) E. Del Nobile,  (2021), 10.1007/978-3-030-95228-0, arXiv:2104.12785 [hep-ph] .
  * Bell  _et al._ (2018b) N. F. Bell, G. Busoni, and I. W. Sanderson, JCAP 08, 017 (2018b), [Erratum: JCAP 01, E01 (2019)], arXiv:1803.01574 [hep-ph] .
  * Abe  _et al._ (2019) T. Abe, M. Fujiwara, and J. Hisano, JHEP 02, 028 (2019), arXiv:1810.01039 [hep-ph] .
  * Ertas and Kahlhoefer (2019) F. Ertas and F. Kahlhoefer, JHEP 06, 052 (2019), arXiv:1902.11070 [hep-ph] .
  * Alexandrou  _et al._ (2020) C. Alexandrou, S. Bacchio, M. Constantinou, J. Finkenrath, K. Hadjiyiannakou, K. Jansen, G. Koutsou, and A. Vaquero Aviles-Casco, Phys. Rev. D 102, 054517 (2020), arXiv:1909.00485 [hep-lat] .
  * Eilers  _et al._ (2019) A.-C. Eilers, D. W. Hogg, H.-W. Rix, and M. K. Ness, The Astrophysical Journal 871, 120 (2019).
  * Bramante  _et al._ (2022) J. Bramante, B. J. Kavanagh, and N. Raj, Phys. Rev. Lett. 128, 231801 (2022), arXiv:2109.04582 [hep-ph] .
  * DeRocco  _et al._ (2022) W. DeRocco, M. Galanis, and R. Lasenby, JCAP 05, 015 (2022), arXiv:2201.05167 [hep-ph] .
  * Page  _et al._ (2004) D. Page, J. M. Lattimer, M. Prakash, and A. W. Steiner, Astrophys. J. Suppl. 155, 623 (2004), arXiv:astro-ph/0403657 [astro-ph] .
  * Yakovlev and Pethick (2004) D. G. Yakovlev and C. J. Pethick, Ann. Rev. Astron. Astrophys. 42, 169 (2004), arXiv:astro-ph/0402143 [astro-ph] .
  * Chatterjee  _et al._ (2022) S. Chatterjee, R. Garani, R. K. Jain, B. Kanodia, M. S. N. Kumar, and S. K. Vempati,  (2022), 10.48550/ARXIV.2205.05048, 2205.05048 .
  * Gardner  _et al._ (2006) J. P. Gardner _et al._ , Space Sci. Rev. 123, 485 (2006), arXiv:astro-ph/0606175 .
  * Neichel  _et al._ (2018) B. Neichel, D. Mouillet, E. Gendron, C. Correia, J. F. Sauvage, and T. Fusco, “Overview of the european extremely large telescope and its instrument suite,” (2018), arXiv:1812.06639 [astro-ph.IM] .
  * Skidmore  _et al._ (2015) W. Skidmore _et al._ (TMT International Science Development Teams & TMT Science Advisory Committee), Res. Astron. Astrophys. 15, 1945 (2015), arXiv:1505.01195 [astro-ph.IM] .
  * Treves  _et al._ (2000) A. Treves, R. Turolla, S. Zane, and M. Colpi, Publ. Astron. Soc. Pac. 112, 297 (2000), arXiv:astro-ph/9911430 .
  * Kouvaris (2008) C. Kouvaris, Phys. Rev. D77, 023006 (2008), arXiv:0708.2362 [astro-ph] .
  * de Lavallaz and Fairbairn (2010) A. de Lavallaz and M. Fairbairn, Phys. Rev. D81, 123521 (2010), arXiv:1004.0629 [astro-ph.GA] .
  * Bertoni  _et al._ (2013) B. Bertoni, A. E. Nelson, and S. Reddy, Phys. Rev. D88, 123505 (2013), arXiv:1309.1721 [hep-ph] .
  * Al-Bakry  _et al._ (2022) M. Al-Bakry _et al._ (SuperCDMS), (2022), arXiv:2203.02594 [hep-ex] .
  * Aprile  _et al._ (2019) E. Aprile _et al._ (XENON), Phys. Rev. Lett. 123, 241803 (2019), arXiv:1907.12771 [hep-ex] .
  * Agnes  _et al._ (2018) P. Agnes _et al._ (DarkSide), Phys. Rev. Lett. 121, 081307 (2018), arXiv:1802.06994 [astro-ph.HE] .
  * Sirunyan  _et al._ (2019a) A. M. Sirunyan _et al._ (CMS), Phys. Rev. D 100, 072002 (2019a), arXiv:1907.03729 [hep-ex] .
  * Aalbers  _et al._ (2022) J. Aalbers _et al._ , (2022), arXiv:2203.02309 [physics.ins-det] .
  * Zhitnitsky (1980) A. R. Zhitnitsky, Sov. J. Nucl. Phys. 31, 260 (1980).
  * Dine  _et al._ (1981) M. Dine, W. Fischler, and M. Srednicki, Phys. Lett. B 104, 199 (1981).
  * Hamaguchi  _et al._ (2019b) K. Hamaguchi, N. Nagata, and K. Yanagi, Phys. Lett. B795, 484 (2019b), arXiv:1905.02991 [hep-ph] .
  * Bramante  _et al._ (2018) J. Bramante, B. Broerman, R. F. Lang, and N. Raj, Phys. Rev. D 98, 083516 (2018), arXiv:1803.08044 [hep-ph] .
  * Bramante  _et al._ (2019a) J. Bramante, B. Broerman, J. Kumar, R. F. Lang, M. Pospelov, and N. Raj, Phys. Rev. D 99, 083010 (2019a), arXiv:1812.09325 [hep-ph] .
  * Bramante  _et al._ (2019b) J. Bramante, J. Kumar, and N. Raj, Phys. Rev. D 100, 123016 (2019b), arXiv:1910.05380 [hep-ph] .
  * McKeen  _et al._ (2021a) D. McKeen, M. Pospelov, and N. Raj, Phys. Rev. D 103, 115002 (2021a), arXiv:2012.09865 [hep-ph] .
  * McKeen  _et al._ (2021b) D. McKeen, M. Pospelov, and N. Raj, Phys. Rev. Lett. 127, 061805 (2021b), arXiv:2105.09951 [hep-ph] .
  * Goldman  _et al._ (2022) I. Goldman, R. N. Mohapatra, S. Nussinov, and Y. Zhang, (2022), arXiv:2203.08473 [hep-ph] .
  * Cortina Gil  _et al._ (2021) E. Cortina Gil _et al._ (NA62), JHEP 06, 093 (2021), arXiv:2103.15389 [hep-ex] .
  * Ahn  _et al._ (2019) J. K. Ahn _et al._ (KOTO), Phys. Rev. Lett. 122, 021802 (2019), arXiv:1810.09655 [hep-ex] .
  * Grygier  _et al._ (2017) J. Grygier _et al._ (Belle), Phys. Rev. D 96, 091101 (2017), [Addendum: Phys.Rev.D 97, 099902 (2018)], arXiv:1702.03224 [hep-ex] .
  * Aguilar-Arevalo  _et al._ (2021) A. Aguilar-Arevalo _et al._ (PIENU), Phys. Rev. D 103, 052006 (2021), arXiv:2101.07381 [hep-ex] .
  * Artamonov  _et al._ (2005) A. V. Artamonov _et al._ (E949), Phys. Lett. B 623, 192 (2005), arXiv:hep-ex/0505069 .
  * Lai  _et al._ (2002) A. Lai _et al._ (NA48), Phys. Lett. B 536, 229 (2002), arXiv:hep-ex/0205010 .
  * Lazzeroni  _et al._ (2014) C. Lazzeroni _et al._ (NA62), Phys. Lett. B 732, 65 (2014), arXiv:1402.4334 [hep-ex] .
  * Abouzaid  _et al._ (2008) E. Abouzaid _et al._ (KTeV), Phys. Rev. D 77, 112004 (2008), arXiv:0805.0031 [hep-ex] .
  * Aaij  _et al._ (2017) R. Aaij _et al._ (LHCb), Phys. Rev. D 95, 071101 (2017), arXiv:1612.07818 [hep-ex] .
  * Bergsma  _et al._ (1985) F. Bergsma _et al._ (CHARM), Phys. Lett. B 157, 458 (1985).
  * Agrawal  _et al._ (2021) P. Agrawal _et al._ , Eur. Phys. J. C 81, 1015 (2021), arXiv:2102.12143 [hep-ph] .
  * Kang  _et al._ (2008) J. Kang, M. A. Luty, and S. Nasri, JHEP 09, 086 (2008), arXiv:hep-ph/0611322 .
  * Jacoby and Nussinov (2007) C. Jacoby and S. Nussinov, (2007), arXiv:0712.2681 [hep-ph] .
  * De Luca  _et al._ (2018) V. De Luca, A. Mitridate, M. Redi, J. Smirnov, and A. Strumia, Phys. Rev. D 97, 115024 (2018), arXiv:1801.01135 [hep-ph] .
  * Gross  _et al._ (2019) C. Gross, A. Mitridate, M. Redi, J. Smirnov, and A. Strumia, Phys. Rev. D 99, 016024 (2019), arXiv:1811.08418 [hep-ph] .
  * Kim  _et al._ (2007) P. C. Kim, E. R. Lee, I. T. Lee, M. L. Perl, V. Halyo, and D. Loomba, Phys. Rev. Lett. 99, 161804 (2007).
  * Perl  _et al._ (2009) M. L. Perl, E. R. Lee, and D. Loomba, Ann. Rev. Nucl. Part. Sci. 59, 47 (2009).
  * Fairbairn  _et al._ (2007) M. Fairbairn, A. C. Kraan, D. A. Milstead, T. Sjostrand, P. Z. Skands, and T. Sloan, Phys. Rept. 438, 1 (2007), arXiv:hep-ph/0611040 .
  * Wells (2003) J. D. Wells, in  _11th International Conference on Supersymmetry and the Unification of Fundamental Interactions_ (2003) arXiv:hep-ph/0306127 .
  * Arkani-Hamed and Dimopoulos (2005) N. Arkani-Hamed and S. Dimopoulos, JHEP 06, 073 (2005), arXiv:hep-th/0405159 .
  * Giudice and Romanino (2004) G. F. Giudice and A. Romanino, Nucl. Phys. B 699, 65 (2004), [Erratum: Nucl.Phys.B 706, 487–487 (2005)], arXiv:hep-ph/0406088 .
  * Preuss and Valencia (2021) C. T. Preuss and G. Valencia, Phys. Rev. D 104, 095030 (2021), arXiv:2109.01261 [hep-ph] .
  * Khachatryan  _et al._ (2016) V. Khachatryan _et al._ (CMS), Phys. Rev. D 94, 112004 (2016), arXiv:1609.08382 [hep-ex] .
  * Aaboud  _et al._ (2019) M. Aaboud _et al._ (ATLAS), Phys. Rev. D 99, 092007 (2019), arXiv:1902.01636 [hep-ex] .
  * Sirunyan  _et al._ (2019b) A. M. Sirunyan _et al._ (CMS), Phys. Rev. Lett. 123, 231803 (2019b), arXiv:1905.10331 [hep-ex] .
  * Sirunyan  _et al._ (2018a) A. M. Sirunyan _et al._ (CMS), JHEP 01, 097 (2018a), arXiv:1710.00159 [hep-ex] .
  * Aad  _et al._ (2021a) G. Aad _et al._ (ATLAS), Phys. Rev. D 103, 112006 (2021a), arXiv:2102.10874 [hep-ex] .
  * Mariotti  _et al._ (2018) A. Mariotti, D. Redigolo, F. Sala, and K. Tobioka, Phys. Lett. B 783, 13 (2018), arXiv:1710.01743 [hep-ph] .
  * Cid Vidal  _et al._ (2019) X. Cid Vidal, A. Mariotti, D. Redigolo, F. Sala, and K. Tobioka, JHEP 01, 113 (2019), [Erratum: JHEP 06, 141 (2020)], arXiv:1810.09452 [hep-ph] .
  * Knapen  _et al._ (2022) S. Knapen, S. Kumar, and D. Redigolo, Phys. Rev. D 105, 115012 (2022), arXiv:2112.07720 [hep-ph] .
  * ATL (2022) (2022).
  * Aaboud  _et al._ (2018) M. Aaboud _et al._ (ATLAS), Phys. Rev. Lett. 121, 211801 (2018), arXiv:1808.02343 [hep-ex] .
  * ATL (2021a) (2021a).
  * CMS (2022) (2022).
  * Misiak and Steinhauser (2017) M. Misiak and M. Steinhauser, Eur. Phys. J. C 77, 201 (2017), arXiv:1702.04571 [hep-ph] .
  * Abe  _et al._ (2020) T. Abe _et al._ (LHC Dark Matter Working Group), Phys. Dark Univ. 27, 100351 (2020), arXiv:1810.09420 [hep-ex] .
  * Aad  _et al._ (2021b) G. Aad _et al._ (ATLAS), JHEP 06, 145 (2021b), arXiv:2102.10076 [hep-ex] .
  * Sirunyan  _et al._ (2018b) A. M. Sirunyan _et al._ (CMS), JHEP 09, 007 (2018b), arXiv:1803.06553 [hep-ex] .
  * Aad  _et al._ (2020a) G. Aad _et al._ (ATLAS), Phys. Rev. Lett. 125, 051801 (2020a), arXiv:2002.12223 [hep-ex] .
  * Bauer  _et al._ (2017b) M. Bauer, U. Haisch, and F. Kahlhoefer, JHEP 05, 138 (2017b), arXiv:1701.07427 [hep-ph] .
  * Haisch  _et al._ (2018) U. Haisch, J. F. Kamenik, A. Malinauskas, and M. Spira, JHEP 03, 178 (2018), arXiv:1802.02156 [hep-ph] .
  * Cepeda  _et al._ (2021) M. Cepeda, S. Gori, V. M. Outschoorn, and J. Shelton, (2021), arXiv:2111.12751 [hep-ph] .
  * Aad  _et al._ (2022) G. Aad _et al._ (ATLAS), Phys. Rev. D 105, 012006 (2022), arXiv:2110.00313 [hep-ex] .
  * Aad  _et al._ (2015) G. Aad _et al._ (ATLAS), Phys. Rev. D 92, 052002 (2015), arXiv:1505.01609 [hep-ex] .
  * Sirunyan  _et al._ (2018c) A. M. Sirunyan _et al._ (CMS), Phys. Lett. B 785, 462 (2018c), arXiv:1805.10191 [hep-ex] .
  * Sirunyan  _et al._ (2019c) A. M. Sirunyan _et al._ (CMS), Phys. Lett. B 796, 131 (2019c), arXiv:1812.00380 [hep-ex] .
  * Aad  _et al._ (2020b) G. Aad _et al._ (ATLAS), Phys. Rev. D 101, 012002 (2020b), arXiv:1909.02845 [hep-ex] .
  * ATL (2021b) (2021b).
  * No (2016) J. M. No, Phys. Rev. D 93, 031701 (2016), arXiv:1509.01110 [hep-ph] .
  * de Florian  _et al._ (2016) D. de Florian _et al._ (LHC Higgs Cross Section Working Group),  2/2017 (2016), 10.23731/CYRM-2017-002, arXiv:1610.07922 [hep-ph] .