Neutron star observations of pseudoscalar-mediated dark matter

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 ~ 0.001 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 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
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) [1,2] have placed remarkable limits on DM over the mass range of m DM ∼ 1 GeV-100 TeV [3][4][5]. The absence of WIMP signals in direct detection experiments has also motivated searches for other DM candidates over a broader range of masses [6][7][8]. 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 DM ∼ 10 MeV in nuclear recoils and m DM ∼ 100 keV in electron recoils [9].
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 [10,11]. 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 DM ∼ 10 −3 [12][13][14][15]. 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 re- * jwcoffey@uvic.ca † mckeen@triumf.ca ‡ dmorri@triumf.ca § nraj@triumf.ca views, see Refs. [16,17]. In particular, capture in neutron stars (NS), beginning with Ref. [18], 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 [19]. 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 , 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 [22,28,36,38,39], DM that scatters primarily inelastically [19,22,24,41], and scenarios where the scattering rate of DM with nuclei is suppressed by positive powers of the DM velocity [21,36,39].
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. [21,36,39]), 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, ultra- 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. violet (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 χ and two realizations of lighter pseudoscalar mediators. The first realization is in the form of an axionlike 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 χ 10 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 χ 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 χ that connects to the SM primarily through a lighter pseudoscalar mediator. We do not attempt to specify the cosmological origin of the χ abundance, but we note that it can be consistent with bounds from indirect detection if its abundance is asymmetric [43][44][45]. 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 [46]. 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 weakscale mediator masses. Thus, taking m a ≤ m χ , we focus on the ALP-mediator benchmark for lighter DM with m χ 10 GeV and the 2HDM+a -mediator benchmark for heavier DM with m χ 10 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 [47][48][49][50]. 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 φ that couples to DM according to This is consistent with a global U (1) P Q symmetry under which [φ] P Q = 1 and [χ L ] P Q − [χ R ] P Q = −1. If φ develops a vacuum expection value (VEV), we can write with the real scalar r(x) and pseudoscalar ALP a(x). This leads to a mass for χ, Changing variables such that χ L → e −iβ L a/f χ L and χ R → e −iβ R a/f χ R with (β L − β R ) = 1 removes the pseudoscalar from Eq. (1) at the cost of introducing new interactions [51]: 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 φ. To be concrete, we assume a set of N Q fermions Q with SM charges (3, 1, 0) together with the coupling [52,53] L ⊃ −y Q φQ R Q L + h.c.
As before, the interaction generates fermion masses m Q = y Q f / √ 2 as well as couplings of the form of Eq. (4) after shifting field variables. The change of field variables also generates an interaction between the pseudoscalar and gluons [54,55]: with c GG = N Q /2. We note that the ALP literature often expresses this coupling in terms of f a = f /2c GG [51]. The interaction of Eq. (6) 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 [54][55][56][57] where f L,R refer to SM fermions and k F and 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 fii −k Fii ). Relative to c GG defined in Eq. (6), we find c ii ∼ 10 −2 c GG . Together with the ALP, this model also contains the radial mode r(x) of the complex scalar, φ(x) = (f / √ 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 ∝ f . The most important couplings of r for our analysis are where the last term is generated by integrating out the heavy color-charged fermions for m Q 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. (6) have been studied recently in Refs. [54][55][56][57], while astrophysical limits from supernova physics are calculated in Refs. [58][59][60]. 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 ∼ f . These are discussed in App. A and can lead to further bounds on the theory, particularly for ALP masses above m a 2 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 10 −2 GeV −1 .

II.2. Heavier Mediator: 2HDM+a
As a benchmark model of a heavier pseudoscalar mediator (m a 10 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 χ [61].

II.2.1. Model and Elementary Couplings
In this model a gauge-singlet pseudoscalar a 0 interacts with the Dirac fermion DM field χ via L χ = y χ a 0χ iγ 5 χ .
The a 0 also connects to the SM by coupling to the H 1 and H 2 doublets of the 2HDM sector according to 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 where V 2HDM is a standard two-Higgs doublet potential [62,63]. On its own, the 2HDM gives rise to physical scalars h and H, a pseudoscalar A 0 , and a charged complex scalar H ± . 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 = v 2 1 + v 2 2 = 246 GeV, tan β = v 2 /v 1 , the CP-even scalar mixing angle α, the 2HDM masses m h , m H , m A0 , m H ± , and the quartic coupling λ 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 where m a and m A are the physical masses of the lighter and heavier pseudoscalars. Together with the angle θ, these three parameters can be used to specify the theory in place of m a0 , m A0 , and µ a12 . After electroweak symmetry breaking, the physical pseudoscalar mass eigenstates a and A couple to DM and SM fermions according to where for a Type-II 2HDM,

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 m A and |θ| 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 β = 3, m H = m H ± = m A = 600 GeV, λ 3 m 2 h /v 2 , and work in the alignment/decoupling limit of cos(β − α) = 0 such that the lighter h scalar couples to the rest of the SM in the same way as the SM Higgs boson [62,63]. 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 [64].
Within this 2HDM-sector benchmark we allow for variations in m a and θ. We briefly summarize the bounds on this parameter subspace here, and discuss them in more detail in Appendix B.
For m a 10 GeV flavor constraints arise. For 1 GeV m a 10 GeV the most important of these come from the loop-induced b-s-a coupling [65,66]. This effective coupling leads to rare decays of bottom mesons, such as B s → µ + µ − and B ± → K ± ( * ) µ + µ − . Limits from the latter process [67] constrain the mixing angle to be smaller than O(few × 10 −3 ) for m a 4 GeV.
When m a ≤ m h /2 the decay channel h → aa can have a significant branching fraction for |θ| 0.05 leading to bounds on θ of this order from searches for exotic Higgs boson decays as well as Higgs rate measurements.
The strongest limits for larger m a ≥ 90 GeV come from LHC searches for H/A → τ τ and imply |θ| 0.1 up to m a 2 m t (for tan β = 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
In this section we provide expressions for the nonrelativistic 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. (6,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. [56,57] (although see Ref. [68] 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 ≡ |q 2 | 500 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 where f π 93 MeV is the pion decay constant. To leading non-trivial order in ε, axion-pion mixing is described by the mixing angle α given by with physical masses ( In terms of these parameters, and working to linear order in ε but to all orders in α (to accommodate m 2 a m 2 π ), we find the physical axion and neutral pion couplings to nucleons N = p, n and DM χ to be with the axial iso-scalar coupling [69,70], the iso-vector coupling g 0 = ∆ p u + ∆ p d 0.339 [70], and where C 0 = c uu + c dd + 2c GG , and ε = ε Θ(m 2 π − m 2 a ).
The couplings in Eq. (20) imply that DM-nucleon interactions are mediated by a combination of axion and neutral pion exchange. The leading Feynman diagrams for spin-dependent N + χ → N + χ scattering are given in Fig. 3(a), and the corresponding matrix element is where with the form factor where σ 3 N = ±1 for N = p, n and the last expression applies for |m 2 a −m 2 π | ε m 2 π . For c uu , c dd → 0, this form factor matches the leading CP-odd gluonic nucleon form factor presented in Refs. [71,72] and given in Eq. (31) below.
In the non-relativistic limit, this matrix element leads to spin-dependent scattering on nuclei with a per-nucleon scattering cross section of where Q 2 = |q 2 | is the 3-momentum transfer, and Q ≤ Q max = 2 µ N χ v χ for incoming non-relativistic DM velocity v χ . A key feature of this result is that the total cross section is suppressed by four powers of v χ after integrating up to Q 2 max .

III.1.2. ALP Spin-Independent
The exchange of the radial mode r leads to nuclear scattering independent of spin and velocity, via the tchannel Feynman diagram in Fig. 3(b). Since the radial mode is heavy, it can be safely integrated out based on the interactions of Eq. (8) to yield an effective operator of the form (9α s /8π) C Gχ χ GG with Matching this gluon operator to nucleons yields an effective SI-nucleon dark matter cross section of where µ N χ is the nucleon-DM reduced mass, 85 GeV is the CP-even gluon nucleon form factor [71,72]. Integrating over Q 2 gives 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 with where and Here MeV, m η = 548 GeV, and the remaining coefficients can be expressed in terms of the quark spin fractions ∆ N q . We use the values from the recent average of 2 + 1 + 1 lattice results of Ref. [70], with a n u,π = a p d,π , a n d,π = a p u,π , a n s,π = a p s,π , and a p,n u,η = a p,n d,η = − 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. [61,65,73] and developed further in Refs. [74,75].
The SI nucleon cross section can be written as with the recent lattice determinations [76] f p u = 0.017 , 4: Contours of constant NS luminosity, logarithmically equispaced, for effective couplings y 2 χ y 2 n = 10 −14 = 1; here Lmax is the maximum NS luminosity that may be imparted by kinetic heating due to incident dark matter flux. Equivalently, these are contours of y 2 χ y 2 n corresponding to NS luminosity Lmax, with values of y 2 χ y 2 n that increase as we go from inner regions of the plot to the outer.
Explicit expressions for the loop-induced C q and C G coefficients in the 2HDM+a are collected in Ref. [74]. 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. (34) 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 β = 3, but at larger values a contribution from two-loop diagrams involving gluons can come to dominate due to a tan 4 β enhancement in loops involving bottom quarks [74]. Higher twist operators are also considered in Ref. [74], 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 [56], and this eliminates the triangle diagram. For the box diagrams, the crucial difference is that the ALP couplings of Eqs. (6,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 UVdependent. The leading contribution to SI scattering in the ALP benchmark is therefore from the radial mode rather than loops.

IV. NEUTRON STAR HEATING
We turn next to study the heating of neutron stars by the infall and capture of local χ 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: giving an escape speed at the surface This configuration is obtained from a Quark Meson Coupling (QMC) equation of state (EOS) of matter near nuclear saturation densities in Ref. [35], 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 →m n , and in our work we set obtained by volume-averaging the NS radius-dependent effective m n shown in Ref. [35]. 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 ρ χ and average DM-NS relative speed v rel (which at the solar position are respectively taken as 0.4 GeV/cm 3 and 350 km/s [77]), the DM mass capture rate is given by [18] is the maximum impact parameter of DM traversing the NS, with 1 + z = (1 − v 2 esc ) −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 where, for optical depth τ , the approximate equality in the first line holds in the optically thin limit considered here. The "capture cross section" above which τ > 1 in the NS core is: where the NS geometric cross section σ 0 The dependence on m χ in Eq. (40) is understood by considering the typical DM recoil energy in the neutron rest frame: 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 ∼Myr [80,81]). Assuming blackbody radiation, the NS luminosity corresponding to a temperature T (in the NS frame) is then This luminosity attains a maximum value L max for unit capture probabilities p σ and p v . Applying these general results to our representative NS configurationwith parameters as in Eq. (35), we find the maximum luminosity L max = 7.6 × 10 24 GeV/s .
This corresponds to a NS temperature seen by a distant observer T = T /(1 + z) of Such a NS temperature is measurable within reasonable integration times at current and imminent infrared 1 As noted in Ref. [78], for mχ 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. [79]. telescope missions [19,21], in particular at the recently launched James Webb Space Telescope (JWST) 2 [83], and the forthcoming Extremely Large Telescope [84] and Thirty Meter Telescope [85].
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 [86]; (b) rotochemical heating in the case of NSs with very small (sub-7 ms) initial spin period for certain nucleon pairing models [27].
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 [87,88] and thereby improve telescope sensitivities [19,21,82]. This effect depends non-trivially on whether or not the DM thermalizes with the NS over the stellar lifetime [31,89], 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 [33,35,39], we generalize our previous results in two ways. First, we extend the nucleon matrix element coefficients C N i of Eqs. (23,25,29,34)by multiplying them by a dipole form factor [71,72], with t = q 2 = −Q 2 and M D 1 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 2 N [71,72]. 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 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 N PS,a → y χ y n /(m 2 a − t) for effective DM and neutron couplings y χ and y n . In Fig. 4 we show contours of constant NS luminosity, corresponding one-to-one to contours of the capture probability p σ , in the plane of m χ versus m a /m χ for the fixed coupling product y 2 χ y 2 n = 10 −14 . Here, the entire region below the contour depicting L max corresponds to p σ = 1. The scaling behaviour of these contours can be understood as follows. For m χ m n , the typical momentum transfer is Q m χ v esc (1 + z) m χ , so we have Combining this with the fact that for m χ < m n Pauliblocking sets p σ ∝ m −1 χ (Eqs. (39) and (40)), we obtain This is indeed the scaling we observe in the m a m χ and m a m χ limits in Fig. 4. For m χ m n , the typical momentum transfer is Q m n v esc (1 + z) m n , and hence In this mass range there is no Pauli-blocking; hence we find these scalings in Fig. 4: We will find the above features reflected in the NS heating sensitivities we estimate below.

V. RESULTS
Having investigated our benchmark pseudoscalarmediated 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 χ . 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 χ /10, representative of a "light mediator" regime, whereas the right panel depicts m a = m χ 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 max given in Eq. (43), where the upper (lower) boundary is for a NS mass of M NS = 1.5 M (2.16 M ). 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. [17]. 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 [35]. Therefore we may infer that, for an observed NS luminosity, there is an uncertainty of O(10) in the minimum ALP-mediated DM cross section deduced from it. In practice we fix this uncertainty to σ uncert = 10, and as the ALP exchange (r exchange) cross section ∝ f −4 ( ∝ f −8 ), the bands spread vertically over factors of σ The solid curves are seen to change slope at two different DM masses: near m χ = 20 MeV in a sharp manner, and near m χ = GeV more gently. The first feature is due to cancellations between the iso-scalar and iso-vector terms in the form factor in Eq. (23) for small momentum transfers. As we raise m χ above m π (while keeping m χ <m n ), Eq. (23) is dominated by the iso-scalar term, and we have σ nχ ∝ |M| 2 /m 2 n ∝ f −4 m 2 χ . In this regime the NS capture cross section σ cap ∝ m −1 χ (Eq. (40)), thus the constraint on 1/f ∝ m −3/4 χ . As we transition to m χ abovem n , we have Q m n , thus the spin-averaged |M| 2 ∝ m 2 χ /f 4 /(m 2 n + m 2 a ) 2 . For m a = m χ , we then have σ nχ ∝ |M| 2 /m 2 χ ∝ f −4 m −4 χ . As the NS capture cross section here is m χ -independent, the constraint on 1/f ∼ m χ . For m a = m χ /10 the turnover occurs at m χ = a few ×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 µ 2 nχ m 2 χ /f 8 dependence on the cross section.
Comparing across the two panels, we also see that near m χ 10 GeV, the constraints on c GG m n /f are about 10 times weaker for m a /m χ = 1 than for m a /m χ = 1/10. This is because in this region the cross section goes as f −4 m −4 a . However, as we go to smaller m χ , and hence smaller m a , the ALP propagator in Eq. (21) 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. (26)). The brown region, which also does not change across the pan-els, is the upper bound from direct detection searches at SuperCDMS [90], Xenon-1T [91], DarkSide-50 [92], and PandaX-4T [5]. 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 O(100 MeV)−O(GeV) due to kinematic limitation in producing mesons. Likewise, for m a GeV a CMS search for the chromo-magnetic dipole moment of the top quark [93] is limited by backgrounds, yielding an upper bound of c GG m n /f 5 × 10 −2 . Thus NS heating could emerge the sole probe in this mass regime for c GG m n /f between ∼ 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 7 × 10 −4 − 3 × 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 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 θ as a function of the DM mass m χ ; we cut off the plot at θ = π/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 χ = 1/3 (m a /m χ = 1). As in the previous sub-section, the bands represent the uncertainty in NS mass-radius configurations, taken to be σ uncert = 10. As the per-nucleon scattering cross section σ nχ ∝ sin 2 2θ (Eq. (27)), the bands spread vertically over a factor of √ σ uncert for small θ.
Further, as σ nχ ∝ m −6 χ for fixed m a /m χ ratios, and as the NS capture cross section is m χ -independent in this regime (see Eq. (40)), the limits on θ ∝ m 3 χ for small θ, 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 [5] and Xenon-1T [4], using the spinindependent loop-induced cross section we derive in Sec. III.2.2. These correspond to exclusion cross sections that range roughly between 4 × 10 −47 cm 2 -4 × 10 −46 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 σ cap 2 × 10 −45 cm 2 ) is weaker than the direct detection limits, thus we do not present it here. We display contours of constant (spin-independent) σ nχ = 10 −47 , 10 −48 , 10 −49 cm 2 , depicting the approximate march of progress expected of direct detection experiments of the next two generations until the so-called neutrino floor is reached [94]. 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 θ are about ten times weaker in the right panel, where a is thrice heavier. This is just as expected from the scaling σ nχ ∝ sin 2 2θ/m 4 a in Eq. (27), implying that the exclusion on θ ∝ m 2 a for small θ. Our NS heating reach is seen to be stronger than the B → K * µµ limits on the mediator, applicable to m a 4 GeV, for m χ /m a = 1/3 in the left panel. In this panel we also outdo collider limits on the mediator from τ τ µµ searches for 12 GeV m χ 25 GeV. Whether we outdo ATLAS limits from undetected final states depends on the unknown mass of the NS observed: a 1.5 M NS does not, whereas greater masses yielding greater NS luminosity (depicted by the red band) do. In the right panel with m χ /m a = 1, our reach is stronger than these limits by only a small margin near m χ = 10 GeV for a 1.5 M NS, while outdoing them for up to m χ 16 GeV as we increase the NS mass to 2.16 M .
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 χ (= m a ) 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 → A → τ τ below m τ τ = 90 GeV. In any case, due to weak dependence on DM mass in the loop-induced SI limits on θ contrasted with the rapid climb (θ ∝ m 3 χ ) in the NS heating reach, we are generally stronger than direct detection for small m χ in the mass range displayed, and weaker for higher m χ . 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 χ 30 GeV for m χ /m a = 1/3 and m χ 15 GeV for m χ /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 KSVZtype 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 [57] and modify the effective nucleon coupling and hence the sensitivity to NS heating. This could be realized in the context of DFSZ-type ALPs [95,96] 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 [25,[36][37][38][39][40].
• 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 [19,22] and other electroweak multiplets [97]), and muon-philic DM in a gauged L µ −L τ model [37]. 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 [7,8,[98][99][100], models that lead to appreciable clustering of DM so that the Earth encounters DM clumps too infrequently [78], 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. [87]), 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 medi-ation (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. [19,21,82] 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 [101][102][103], 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. 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. [57] 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 → T 2 (r) = 1/2. The key phenomenological divider here is the QCD scale ∼ m π . A non-zero c GG induces sizeable ALP couplings to the photon via the QCD anomaly, c γγ −1.92c GG [55], so that for m a < m π the diphoton mode a → γγ dominates the ALP decay width. For m a > m π the branching to hadrons dominates, but decays to muons a → + − at the sub-percent level turn out to be important.
For m a < m π the strongest terrestrial constraints come from searches for K + → π + + X inv at NA62 [104], with X inv an invisible state. In this range of m a the a → γγ decay length exceeds the 10 m detector size and is thus unseen [57]. Limits from K L → π 0 X inv at KOTO [105] are weaker due to suppressed CP-conserving couplings, and even weaker limits are set by B → K * X inv at Belle [106] and π + → ae + ν e at PIENU [107]. For somewhat large m a and c GG the a → γγ mode becomes visible, and limits are obtained from K + → π + γγ and K L → π 0 γγ at E949, NA48, NA62, and KTeV [108][109][110][111].
For 2m µ < m a < m B limits from B + → K + a(µµ) and B 0 → K * a(µµ) at LHCb [67,112] apply, and for m a GeV the strongest limits come from those placed by CMS [93] on the chromo-magnetic dipole moment of the topμ t , defined by the effective operator The ALP contributes toμ 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 ∝ c 2 GG [57]. For Λ QCD m a < 800 MeV we have limits from the visible mode a → γγ at the CHARM beam dump [113], recast by F. Kling in Ref. [114]. Near m a = m η 0 500 MeV, due to strong ALP-η 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. [58] 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, γp → ap, however the limits from this process assuming non-zero c γγ (trivializing all other couplings) turn out to be much weaker than those from nucleon bremsstrahlung [58], 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" [60]. 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. [59] by considering the boost in explosion energy from a → γγ in several low-energy supernovae observed. We note that other astrophysical limits apply for coupling ranges below our description above in the m a 100 MeV region [60], 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 (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 10 TeV [115][116][117][118]. After QCD confinement, the heavy quarks form bound colorneutral states with themselves as well as with ordinary quarks. The analysis of Ref. [117] argues that rearrangement reactions convert the vast majority of Q-SM mixed states to more deeply bound, neutral QQQ and QQ 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-SM /n B 10 −21 from searches for fractional charges in oil drop experiments using standard and meteoritic materials [119,120]. 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 Qs was created.
Direct collider searches for stable Q quarks place bounds that are independent of cosmology. Following QQ pair production, the heavy quarks would be expected to hadronize and lose energy as they pass through the detector [121], in analogy to long-lived R-hadrons considered in the context of (very) split supersymmetry [122][123][124]. Generalizing the calculations of Ref. [125] to colortriplet fermions, we estimate that the LHC searches of Refs. [126,127] translate into m Q 1500 GeV. Note that since m Q /(f /c GG ) = y Q N Q /2 √ 2, this can be achieved for c GG /f 10 −2 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 ≤ 10 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 GeV, the leptophobic Z search of Ref. [128] translates into the bound c GG /f 3 × 10 −2 GeV −1 , but we do not know of any useful bounds at lower masses. The heavier radial mode with m r ∼ f decays primarily via r → 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. [129] gives no bound for N Q = 1 and c GG /f 2 × 10 −2 GeV −1 for N Q = 10 and m r = f . We estimate that monojet-like searches for r → χχ decays such as Ref. [130] 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 For the ALP, this corresponds to g aγγ = 2α c GG /(3πf ) 1.7 × 10 −3 c GG /f in the standard notation [51]. 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 [56]. When both photon and gluon interactions are present, as we have here, bounds that rely on on-shell a → γγ are weakened further by competing decays to hadrons while collider limits can become much stronger [131][132][133]. Extrapolating the recent LHC boosted diphoton search of Ref.
[134] to this scenario, we find c GG /f 5 × 10 −4 GeV −1 for m a = 10 GeV. Unfortunately, this study does not extend to lower ALP masses and therefore we do not obtain a related limit in the m a ∈ [2, 10) GeV mass range of primary interest. The data scouting method proposed in Ref. [133] is projected to have sensitivity to c GG /f ∼ 2 × 10 −3 GeV −1 . We have also studied the impact of the induced couplings of the radial mode to photons of Eq. (A2), but we do not find any limits that are stronger than the dijet bounds discussed already because of the dominance of r → aa decay mode together with the very small branching ratio for a → γγ. Finally, we note that direct searches for pair-produced heavy quarks give bounds that are similar to or slightly weaker than m Q 1500 GeV [135][136][137]. mostly singlet. As a benchmark for our study, we fix m A = m H = m H + = 600 GeV, tan β = 3, λ 3 m 2 h /v, and cos(β − α) = 0 corresponding to the so-called alignment/decoupling limit [62,63]. The dark matter dynamics then depend on m a and θ, 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 β = 3, the heavy Higgs masses are large enough to be consistent with flavour bounds from B → X s γ [138] and B s → µµ [74,139], as well as collider searches such as H + → tb [140] and H/A → ττ [141,142]. 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 [143]. Similarly, the alignment limit yields an SM-like Higgs boson with production and decay rates effectively identical to the SM [62,63].
While the 2HDM sector is largely safe in this benchmark, the presence of a lighter pseudoscalar implies bounds on the mixing angle |θ| 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 10 GeV. In our region of parameter space, these are dominated by the flavorchanging coupling of a to a bottom and strange quark induced at one-loop [65]. Using the results of [66], we compute the rate for B ± → K ± ( * ) (a → µ + µ − ) and compare to the LHCb search for promptly-decaying bosons in this channel [67]. In addition, there is a subdominant constraint that comes from requiring that the branching fraction for B s → µ + µ − agree with the world average of (3.01 ± 0.35) × 10 −9 [69]. The resulting limits on θ from these processes are shown in Fig. 6.
Pseudoscalar masses m a m h /2 are also strongly constrained by searches for rare SM Higgs decays. The decay width for h → aa in our benchmark is [61,74,143] With m A = 600 GeV this width exceeds the total SM Higgs width for |θ| 0.05 and thus Higgs measurements imply bounds on the mixing angle near this level. Since we focus on m a ≤ m χ , 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. [143,144]. Exotic decays of the SM Higgs have been searched for the by LHC collaborations [145], and in Fig. 6 we show estimated limits from ATLAS searches for h → aa → bbµµ [146] and h → aa → τ τ µµ [147] as well as CMS searches for h → aa → bbτ τ [148] and h → aa → 4µ [149]. 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 und < 0.15 [150,151]. We also show the corresponding limit on |θ| 0.04 in Fig. 6 under the assumptions of unmodified SM Higgs production and that h → 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 χ with dominant a → χχ decays [65,139,144,152], 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 → τ τ can be reinterpreted to give the strongest limits on this scenario. For our benchmark of tan β = 3, the most important production channel is in association with bottom quarks. We estimate the production cross section by rescaling the bbφ cross sections tabulated in Ref. [153] by the appropriate factor of (tan β sin θ) 2 , and derive bounds on |θ| by comparing to the limits obtained from φ → τ τ searches at ATLAS [142] and CMS [141]. The exclusions we find based on these studies are shown in Fig. 6. by [30] p v = where we set v halo = 10 −3 c and v 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 CPeven 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 O(1) fraction of scattering events to evade capture for m φ /(2µ nχ ) ≤ 0.01, as also shown in Ref. [30].