Probing Light Dark Matter with a Hadrophilic Scalar Mediator

We investigate the thermal cosmology and terrestrial and astrophysical phenomenology of a sub-GeV hadrophilic dark sector. The specific construction explored in this work features a Dirac fermion dark matter candidate interacting with a light scalar mediator that dominantly couples to the up-quark. The correct freeze-out relic abundance may be achieved via dark matter annihilation directly to hadrons or through secluded annihilation to scalar mediators. A rich and distinctive phenomenology is present in this scenario, with probes arising from precision meson decays, proton beam dump experiments, colliders, direct detection experiments, supernovae, and nucleosynthesis. In the future, experiments such as NA62, REDTOP, SHiP, SBND, and NEWS-G will be able to explore a significant portion of the cosmologically motivated parameter space.


I. INTRODUCTION
The dark matter puzzle is among the most pressing problems in particle physics and cosmology today. In recent years, a growing experimental and observational program to search for non-gravitational dark matter interactions has coincided with a broader theoretical exploration of viable dark matter models, novel cosmological histories, and associated phenomenological signatures. Through these efforts, the paradigm of a light dark sector containing new (sub-)GeV-scale singlet dark particles coupled to ordinary matter through a light mediator has emerged as a compelling possibility and remains under active investigation.
An excellent survey can be found in the recent community studies [1,2].
The experimental signatures of a particular dark sector model are, to a large extent, governed by the couplings of the mediator particle to ordinary matter. For instance, in the popular vector portal dark matter model [3][4][5][6], the "dark photon" mediator couples to electrically charged particles via kinetic mixing with the ordinary photon [7,8]. A rich phenomenology results due to the couplings to both leptons (e.g., electron beam fixedtarget experiments, low-energy e + e − colliders, direct detection via electron scattering, lepton anomalous magnetic moments) and quarks (e.g., proton beam fixed-target experiments, direct detection via nuclear scattering, rare meson decays); see Refs. [1, 2] and references therein. Besides kinetically mixed dark photons, a variety of mediators have been explored in the literature, including Higgs portal scalars [9,10], neutrino portal fermions [11][12][13][14][15][16][17], and vector bosons coupled to anomaly free currents [18][19][20][21][22][23], all of which have distinctive patterns of couplings to SM particles and resulting phenomenology. Such wide theoretical investigation is required to identify the full range of phenomena associated with dark sectors and discern the physics potential of proposed searches and new experiments with respect to existing constraints over a broad range of models and couplings.
It is in this context that we are motivated to explore dark sectors with hadrophilic (or leptophobic) mediators. Finding viable constructions with light hadrophilic mediators is challenging in comparison to the mediators mentioned above. Leptophobic vector mediators coupled to, e.g., baryon number, have been considered in the past [24][25][26][27][28][29][30][31] and are possible in principle. However, it was recently shown that such scenarios face stringent constraints due to enhanced production of the longitudinal mode in a variety of rare decays [32,33], a result that can be traced to the anomalous nature of the baryon number symmetry. Lep-tophobic scalar mediators, on the other hand, face a different set of challenges, including new flavor-changing neutral currents (FCNC) and a naturalness problem related to the light scalar mass. Typically, a nontrivial flavor hypothesis (e.g., Minimal Flavor Violation [34] or alignment) on the scalar-quark couplings is needed to satisfy FCNC constraints, while naturalness arguments would point to weaker couplings for smaller scalar masses. An extensive analysis of these issues was recently carried out in Ref. [35], where it was found that "flavorspecific" scalar couplings, i.e., scalars that couple to a specific quark mass eigenstate, can satisfy existing FCNC constraints even for relatively sizable couplings in the natural regions of parameter space. This framework therefore provides a promising point of departure to study a light leptophobic dark sector, and this is the task undertaken in this paper. We note that a similar analysis of the flavor-alignment hypothesis in BSM theories is presented in Ref. [36]. In addition, previous work considered leptophobic scalars that couple dominantly to gluons or to top quarks in the sub-GeV range [37]. Here we will instead focus on the case in which the scalar couples dominantly to the up-quark, finding some differences in comparison to the gluon-specific and top-specific scenarios that will be highlighted below.
In particular, we find that a diverse and complementary set of experimental and observational probes, including precision meson decay measurements, beam dump experiments, colliders, direct detection, supernovae, and nucleosynthesis, already provides significant coverage of this scenario. Nevertheless, there are viable regions of parameter space where thermal freezeout can set the correct relic dark matter abundance, both through direct annihilation to hadrons as well as through secluded annihilation to scalar mediators. We also describe several promising ongoing or future experiments that will be able to further test this scenario, including NA62, REDTOP, SHiP, SBND, and NEWS-G.
The remainder of this paper is organized as follows. In the next section, we define the basic framework for the dark sector with up-quark specific scalar mediator couplings. In Section III, we investigate the phenomenological implications of this scenario for meson decays, beam dump experiments, colliders, direct detection, astrophysics, and cosmology.
Section IV contains our conclusions. We have also included an Appendix which presents a description of the hadronic couplings for a general flavor-diagonal scalar mediator.

II. FRAMEWORK
The theoretical framework for a flavor-specific scalar mediator has been presented in Ref. [35], and we begin by reviewing its essential features. A singlet scalar mediator S with mass m S is assumed to couple the visible and dark sectors. We will take dark matter to be a SM singlet Dirac fermion χ with mass m χ that is charged under a Z 2 stabilizing symmetry. The mediator is assumed to couple dominantly to up quarks through a dimensionfive operator generated at a UV scale M . The relevant terms in the Lagrangian are where in the second line we have defined the effective coupling with v = 246 GeV the Higgs vacuum expectation value (vev). Other possible dimension-five operators involving S,Q L ,U R , such as ∂ µ SŪ R γ µ U R and iSŪ R / DU R + h.c. , are related to the one written in Eq. (1) by a field redefinition. Besides these, additional operators such as SG µν G µν are expected on general grounds, but their presence and size depends on the particular UV completion. We will assume that the c S coupling in Eq. (1) dominates over other dimension-five operators.
viewed as a spurion, transforms as c S ∼ (3,3, 1). The assumption that S couples dominantly to the physical up quark, such that in the mass basis c S ∝ diag(1, 0, 0), implies that c S breaks the flavor symmetry according to the U (1) u factor is aligned with and identical to the one left unbroken by the up-quark Yukawa spurion. We emphasize that this symmetry breaking pattern is a hypothesis on the form of our low energy effective theory. It would certainly be worthwhile to explore UV model constructions which realize alignment and flavor-specific structures; for some promising directions along these lines, see Refs. [38,39]. Setting aside questions on the UV origin of such structures, the vast majority of dangerous FCNC processes are suppressed to acceptable levels once this flavor hypothesis is made. However, as we will see in Section III, rare flavor changing meson decays (e.g., kaon decays) can still provide a sensitive probe if the scalars are light.
We now turn to a discussion of the scalar potential for S. As is well-known, a light scalar with large couplings to heavy UV states is expected to receive large radiative contributions to its mass (quadratic) and, in the case of a singlet scalar, tadpole (linear) terms in the potential. A standard naturalness argument would then suggest that the lighter the scalar particle is, the smaller its couplings should be. This argument applies without any ambiguity to c S in Eq. (1); since it is a dimension-five operator we must introduce degrees of freedom at the UV scale M , which would give a large correction to the scalar mass unless c S is sufficiently small. In addition, in the case of the singlet considered here, the tadpole induces a vev for S and in turn a contribution to the effective up-quark Yukawa coupling. As discussed in detail in Ref. [35], the size of the leading two-loop corrections to the scalar mass and up-quark Yukawa (as well as other terms in the potential) are small provided We note that the criterion in Eq. (3) is derived considering only the interactions and fields in the effective field theory (1). In particular, it often happens that in explicit UV completions there is a larger one-loop correction to the scalar mass which can result in a moderately more restrictive criterion than the one given in Eq. (3). In any case, it is worth emphasizing that one can of course explore regions of parameter space not satisfying (3) at the expense of fine-tuning. For this work, we will use (3) to provide a rough picture of the boundary between the natural and tuned regions of parameter space.

A. Thermal Cosmology
For m χ < m S , and assuming there are no lighter states in the dark sector, χ will annihilate directly to SM particles through an s-channel mediator. In this case, there is a region in the (m χ , g χ , g u ) parameter space that predicts the observed dark matter relic density assuming a standard thermal cosmology. We take this region as a motivated target in parameter space, though we note that the χ density may also be set by non-thermal mechanisms, such as an early χ-χ asymmetry. If g χ has a non-zero phase, s-wave annihilation can occur, which is constrained for symmetric dark matter at the thermal level by Planck measurements of the Cosmic Microwave Background (CMB) for m χ < ∼ 10 GeV [40] and by Fermi-LAT dwarf spheroidal galaxy observations for m χ < ∼ 100 GeV [41]. As we will be interested in m χ < ∼ O(GeV), we henceforth assume that g χ is purely real, in which case the cross section for annihilation to SM final states is given by For real g χ we observe that the annihilation rate is velocity suppressed, and the CMB and gamma-ray constraints are trivially satisfied. We note that Γ S in Eq. (4) is the full S width while Γ S m S =2mχ is the width of the scalar evaluated at m S = 2m χ . For m S well above Λ QCD , we can use the perturbative width for S → uū, which is given by Γ S→uū = 3g 2 u m S /8π. At lower scalar masses, we must take hadronic effects into account to calculate Γ S ; our procedure will be described below while the full details are given in the Appendix. Finally, for the thermally averaged annihilation cross section, σv rel , we use v 2 rel = 6/x, with x = m χ /T and T the temperature.
In the opposite regime, m χ > m S , dark matter will dominantly annihilate to pairs of mediators ("secluded" annihilation [4]) for most parameter choices. For purely real g χ , and in the limit m S m χ , the annihilation cross section is given by We note that this cross section is also velocity suppressed and thus safe from CMB and gamma-ray constraints. It also only depends on the dark matter -mediator coupling g χ , and the correct thermal abundance can be achieved provided the mediator has only a minuscule coupling g u with the SM to maintain kinetic equilibrium with the bath.
With real g χ as motivated above, even a small imaginary component of g u will generate a large neutron EDM that is generically in tension with the experimental limit [42]; see Refs. [35,43] for detailed discussions. We will therefore take g u to be real as well from here on so that the scalar interactions respect CP . While thermal DM annihilating through a real scalar interaction generally faces strong spin-independent direct detection bounds, we will be interested in the sub-GeV regime where these constraints are relatively weak. A full discussion of this issue will be postponed to Section III.

B. Hadronic couplings of S
At scales between Λ QCD and the cut-off scale M of the dimension-five operator in Eq. (1), the above Lagrangian provides a good description of the theory. Below the QCD scale, however, the interactions of S are naturally written in terms of its hadronic couplings. Many of the probes of the sub-GeV dark sector involve sub-GeV momentum transfer, so it is important to establish such a description. This is analogous to studies carried out long ago for a light SM Higgs boson [44][45][46][47][48][49] and more recent studies of light Higgs portal scalars (see, e.g. [50][51][52]). We will begin by discussing the use of the chiral Lagrangian to estimate the where Σ = e 2iπ/f contains the pNGB fields, π = π a T a , f ≈ 93 MeV is the pion decay constant, and the scalar couplings of the quarks are contained in The spurion χ contains both the usual quark masses as well as the effective coupling of S to quarks. In the above, B is a dimensionful parameter that can be determined by expanding the Lagrangian to obtain the physically observed meson masses, yielding B GeV. The Lagrangian (6) contains the meson mass terms and trilinear interactions of S with the mesons 1 . As expected, the Sππ couplings are all of order g u B.
We have also included the η meson in our description; see the Appendix for details.
At higher energies, chiral perturbation theory is insufficient to describe the interactions between S and the mesons. We parametrize the couplings by form factors [48,50,51], e.g.
for the pions π(p)π(p )|m uū u + m dd d|0 = Γ π (s), as a function of the momentum transfer squared s = (p + p ) 2 . In practice, the form factor We are now ready to consider the decays of the mediator S, which are needed as input to the phenomenological analysis presented in Section III as well as the dark matter annihilation cross section (4). For m χ < m S /2, the decay S → χχ occurs, and without a strong bound on g χ we simply assume that it dominates the decays of S. On the other hand, for m χ > m S /2, S decays exclusively to photons or hadrons, depending on its mass relative to the pion threshold m S = 2m π . For a scalar light enough to have no available tree level decay modes, the only available decay is to photons. The partial width for this decay is where the loop function is The sum runs over all possible quarks in the loop, and when a first generation quark dominates, we use the constituent quark massm u =m d = 350 MeV.
For m S above the pion threshold, S decays to hadrons. As for the case of DM annihilation, we must use the proper form factors for the interactions of S with the mesons when m S is near Λ QCD . Again following [51], we match the partonic decay width for S → uū at high m S to the sum of S decays to mesons at low masses. Further details are contained in Appendix A. Fig. 1 shows the resulting branching ratios and decay lengths of the scalar.
For m S < 2m π , the only possible decay of S is through a loop to photons. Above the pion threshold, S → ππ is the only significant decay until kaons are kinematically accessible, strongly affecting S decay through the f 0 (980) resonance. For m S > ∼ 1.3 GeV, we start to lose predictive power on the individual hadronic decay channels of S as the "other hadronic" decay channels, which we do not compute explicitly, are dominant. Nevertheless, we expect that our estimate for the S total width is still reasonable. This calculation will have an impact in the next section, where the lifetime of S is important in determining the relevant range of couplings g u that may be probed by a given experiment for fixed m S .

III. PHENOMENOLOGY
In this section, we explore the viable space for MeV-GeV scale dark matter with an upphilic scalar mediator. We describe current limits and future prospects from meson decays, beam dump experiments, Big Bang Nucleosynthesis (BBN), supernovae, direct detection, and colliders.
We will consider both the cases where the scalar decays visibly to SM particles, as described in Section II, as well as where the scalar decays invisibly to χχ. We first note that in the former case, there is a sharp change in the lifetime of S at m S = 2m π as indicated by the right panel of Fig. 1, due to the only available decay channel for lighter S being the loop level decay to γγ. For m S below the pion threshold, it is clearly possible that the scalar decays to SM states, but lives long enough that this decay is not observed in laboratory experiments, and so the scalar effectively appears to be invisible. Whether the photons from S decay are visible in a given search depends on the detector setup as well as the energy with which the scalar is produced. Generally, these photons can be seen either if the coupling is large such that the loop decay proceeds rapidly, or if the scalar is produced with some boost and the detector is sufficiently far away from the point of production. Below, we will show examples of both types of searches which can directly observe the decay S → γγ. Otherwise, searches that would normally apply to an invisibly decaying S can still be used. Here, the requirements are the converse: the coupling must be low enough or the detector must be sufficiently small that S decays outside of it. Above the pion threshold, the decay of S to pions generally proceeds promptly except at very low couplings, where S can still be longlived. When the decay is prompt, the S can in principle be detected as a pion resonance.
In practice, we will see that there is only a small region of m S where these decays provide a useful constraint on g u .
We will also discuss the scenario where the new scalar decays invisibly, assuming that the branching fraction of S → χχ is approximately 1, and furthermore that this decay occurs promptly. This is expected if m χ < m S /2 and the coupling g χ is larger than the SM coupling of S. In the following, when considering a possible invisible decay mode of the scalar, we will assume that it dominates.
In principle the couplings g u , g χ and the masses m χ , m S are free parameters in our scenario. However, beyond the distinction of whether S decays visibly or invisibly, the mass ratio m χ /m S has little practical impact on many searches for S that involve only the Sūu coupling, such as those at fixed-target experiments. As long as g χ is sufficiently larger than g u , changing it does not lead to significantly different phenomenology 3 . With this in mind, we will generally show limits in the m S -g u plane. We note that for m χ < m S the thermal DM target and direct detection limits do not occupy a unique position on this plane, owing to the freedom to choose m χ and g χ . On the other hand, for m χ > m S , the thermal annihilation scenario is completely independent of the SM couplings of S since secluded annihilation χχ → SS depletes the χ abundance; thus the correct relic density can be obtained anywhere in the m S -g u plane, with an implied relation between m χ and g χ from Eq. (5). With these points in mind, we now explore the phenomenology of our scenario, with an eye towards both differences from the usual Higgs-like scalar case and opportunities for improvement in the future.

A. Meson decays
For sub-GeV m S , there are many possible meson decays that contain S in the final state, and through these decays S can be copiously produced at precision decay experiments.
Unlike a Higgs-like scalar, for which meson-induced scalar production typically involves the top quark coupling in a penguin diagram, a scalar which couples predominantly to the first generation is produced directly from the quark content of the mesons. We estimate the production of scalars from the decays of η and K mesons using the chiral Lagrangian of Eq. (6). 3 In the region 2m χ ≈ m S , the relic density calculation is affected significantly by resonant annihilation.
The resulting branching ratios for the leading decays of the mesons which produce S are We briefly comment on the production of S from heavy meson decays, e.g., B → KS, where the chiral Lagrangian approach is not applicable. This is typically one of the largest production channels for a Higgs-like scalar, but turns out to be irrelevant for an up-philic scalar. We have estimated two possible contributions to this channel. First, we have computed the tree level weak decay with emission of S from the spectator up-quark, dressed with meson form factors and wave functions parametrizing the momentum distribution of the quarks within the mesons [53]. However, due to the significant V ub suppression in the amplitude, we find a branching ratio that is 1-2 orders of magnitude below current limits [54] even for order one g u . Second, we note that in principle S has an induced interaction with gluons at two loops, which in turn feeds into an Stt coupling, giving penguin diagram contributions to the decay. However, the combination of the S shift and up quark chiral symmetries implies that any SG µν G µν interaction be proportional to both g u and Y u . Together with the loop suppression, the resulting contribution to B decays to the scalar is far too small to be of any importance. This behavior is in contrast to scalars whose interactions are dominated by top quark and gluon couplings, where searches for rare heavy meson decays set relevant bounds [37]. Given these estimates, we do not consider limits from B decays further.
We now turn to specific experiments which can exploit the meson decays of Eq. (11). In many cases, there are direct measurements of exclusive meson decay channels. First, the decay η → πγγ has been measured at the Mainz Microtron (MAMI) with 10% precision [55].
We use this measurement to constrain the decay η → πS, S → γγ [56], conservatively requiring only that the branching ratio of such decays which occur within the Crystal Ball detector be smaller than the 2σ upper bound on BR(η → πγγ) of 3.0 × 10 −4 , i.e., making no assumption on the SM contribution to this final state. To calculate the fraction of the decays which produce photons within the detector, we use the fact that the η is produced nearly at rest within the detector of radius 25 cm [57], which in our scalar mass region of interest puts a lower bound on g u .
In addition, 4.7×10 6 η → π 0 π + π − decays have been recorded by KLOE [58], with a Dalitz analysis that shows no signs of resonances as would be produced if there were a contribution from η → π 0 S, S → π + π − . For a given value of m S , we estimate the background in the corresponding m S = m π + π − bin by fitting the data in the other bins using a quadratic polynomial in m π + π − to describe the SM η → 3π matrix element. We then require that the background plus the contribution from η → π 0 S not exceed the 3σ upper limit in this bin. For the values of g u probed here, the S is very narrow and does not populate multiple m π + π − bins. This limit is valid in the mass range 2m π < m S < m η − m π . The proposed REDTOP experiment [59] hopes to produce around 10 13 η/year using a 1.9 GeV p beam on a beryllium target. Given the branching ratio for η → π 0 π + π − of 23% this could result in a sample of π 0 π + π − events about 10 5 times larger than at KLOE, which could allow the limit on g u to be improved by a factor of ∼ 20.
For an invisibly decaying scalar, the analogous η decay channel η → πS, S → χχ could be considered. However, there is currently no search for η → π + invisible. Performing such a search at the upcoming REDTOP experiment would be challenging since the signature would be simply two photons (that reconstruct a π 0 ) recoiling against missing momentum.
However, because the longitudinal momentum of the S is unknown, reconstructing the η would be challenging. A second phase at REDTOP is also envisioned, using a higher energy beam to produce around 10 11 η /year. The branching for η → π 0 π 0 η is 23% and the extra kinematical constraints from the π 0 → γγ decays could allow the "tagging" of an η that decays partially invisibly. A branching of Br(η → π 0 S) 10 −7 , corresponding to g u 10 −6 , would lead to O(10 3 ) such events. A detailed study of the prospects at REDTOP for this decay mode would be highly desirable.
Precision kaon measurements can also constrain our scalar, through the decay K → πS.
In particular, experiments E787 and E949 at Brookhaven searched [60][61][62][63] for the decay K → πνν, which constrains the rare decay K → πS where S is either long-lived or decays invisibly. When S decays visibly but m S < 2m π and g u is sufficiently small, S does not decay within the E949 detector, and we calculate the expected number of scalars from K → πS which would produce such "invisible" signatures by asking how many scalars would decay to photons outside the 1.45 m radius of the E949 detector, as a function of g u and m S .
The upper edge of the exclusion region in the m Sg u plane in the left panel of Fig. 2 corresponds to the coupling where the S decays inside the detector, so that the scalar is no longer effectively invisible. For a genuinely invisibly decaying scalar, of course, the same search applies without an upper limit on g u .
The K → πνν mode will be measured to 10% precision by the NA62 experiment [64], which uses the 400 GeV CERN SPS proton beam to produce a secondary 75 GeV kaon beam, allowing the study of order 10 13 decaying kaons over the life of the experiment. To estimate the increased sensitivity, we scale the E787/E949 result by the increase in the precision on the SM branching ratio for the decay being studied, following similar projections for axionlike particles [65]. The actual eventual reach of NA62 will also depend on the theoretical uncertainty in the SM prediction for the decay branching ratio, which is currently 10% [66]; we neglect this here.
Lastly, the decay η → πS can be used to constrain even heavier scalars than those probed by η or kaon decays. Due to the strength of the other meson bounds, we focus only on the mass region m K − m π < m S < m η − m π , above the kinematic reach of kaon decays. If it decays to SM particles, a scalar of this mass would contribute to the decay η → 3π, which has been measured by BESIII [67]. We require that for such a scalar, the branching ratio for η → πS not exceed the measured branching ratio for η → 3π. On the other hand, if S decays invisibly, we impose that the partial width Γ(η → πS) not exceed the PDG fit for the total η width [68] of 196 keV. Additionally, a high energy REDTOP run could increase the number of reconstructed η substantially.

B. Beam dumps
Having discussed precision decay measurements, we now describe the impact of beam dump experiments that can search for S. At proton beam dumps, because of the significant production cross section for η mesons and its narrow width, η → πS can provide a significant scalar production channel. Furthermore, the η production rate at proton beam dumps is typically less than that of η by about an order of magnitude, but allows for additional access to the scalar mass region m η − m π < m S < m η − m π ; note that as this region is above the pion threshold, the additional gain is at very low coupling. Although m K and m η are similar, the branching fraction of kaons to scalars is smaller than that of η to scalars at the same coupling, and kaons are produced less copiously in proton beam dumps. Therefore, kaon-induced production of S is comparatively small, and we neglect it below. We emphasize again the difference with the case of a Higgs-like scalar: there, B and K mesons are the usual meson decay sources of production at beam dumps. However, for a scalar that preferentially couples to the first generation, the effective coupling of the heavy mesons to the scalar is very small, as discussed above, and the primary sources of scalars from meson decays are instead the light mesons. At proton fixed-target experiments, hadrophilic scalars can also be produced through bremsstrahlung, but we consider only meson decays here. It would be useful in the future to perform a detailed calculation of the additional bremsstrahlung production mode, which would improve the limits shown here.
The CHARM collaboration performed a search for axion-like particles produced from the 400 GeV SPS proton beam, decaying to γγ or leptons producing electromagnetic showers, in a detector located 480 m downstream and approximately 10 mrad off axis [69]. We recast this search to place a limit on scalars produced from η decays. To obtain the total S production rate, we use CHARM's estimate of pion production, and then apply a scaling factor obtained from a simulation in CRMC [70] showing that approximately one η is produced for every ten neutral pions within the geometrical acceptance of the CHARM detector 4 . The distribution of the η energy affects the boost of the scalar produced in η → πS, and we simply assume that all S particles are produced with energy 25 GeV, the average of the E η spectrum [69], as in previous studies [71]. Combined with Br(η → π 0 S) above and the lifetime estimates for S from the previous section, we may then calculate the number of γγ events that would have been seen by CHARM for m S < 2m π . As no events were seen, we simply require that no more than three γγ events would have occurred in the CHARM detector, assuming perfect reconstruction efficiency. Above the pion threshold, we note that each decay S → π 0 π 0 would produce 4γ events, which would have again been visible at CHARM as electromagnetic showers. We can thus still place a bound from the CHARM search for m π > 2m S , occupying a region at much lower coupling due to the tree level decay of the scalar. Future proton beam dumps will also be sensitive to S production from mesons. We consider the case of SHiP [71] in complete analogy with CHARM. SHiP would again use the 400 GeV SPS proton beam to produce very weakly coupled light particles, though with a much closer detector than CHARM, only about 70 m from the interaction point. Again we use CRMC to simulate η( ) production from proton interactions with the target material, and employ the same assumptions on the scalar energy spectrum and reconstruction efficiency as for CHARM above. Unlike for CHARM, the η production is sufficient to probe an additional part of parameter space in the relevant mass regime at g u ∼ 10 −8 . It would also be interesting to examine the sensitivity of other proposed detectors at the LHC and beyond targeting long-lived particles, including CODEX-b [72], FASER [73], MATHUSLA [74], and SeaQuest [75,76].
While thus far we have only discussed proton beam dumps, the loop-level coupling of our scalar to photons implies that S may be produced through Primakoff production at electron beam dumps as well, analogously to the muon-philic case [35]. We have checked the limits on the loop-level interaction from searches for axion-like particles at experiment E137 at SLAC [77][78][79]. The resulting bound largely overlaps with that from CHARM, and does not constrain any additional parameter space once the other limits in this section are taken into account. Now, when S decays to dark matter, the standard beam dump tests searching for the decay products of new particles that are produced in fixed-target collisions and decay far downstream are no longer applicable. However, through the decay of S, light DM can be produced in the primary collisions of protons in a beam dump and subsequently detected through its scattering with nucleons in a downstream detector [25,27,28,30,31,[80][81][82][83][84][85][86].
This approach has been employed recently by the MiniBooNE-DM collaboration [87][88][89], and the null result from their search in the nucleon elastic scattering channel [87] leads to relevant constraints on our scenario. To estimate the potential dark matter scattering yield and derive the constraints implied by the search [87] we have performed a Monte Carlo simulation using methods similar to those employed in the BdNMC package [85]. During the dedicated MiniBooNE-DM run, 1.86 × 10 20 protons-on-target (POT) from the Fermilab Booster (8 GeV kinetic energy) were directed onto an iron beam dump. Focusing on the regime m S > 2m χ , we consider the following DM production chain: 1) Meson production in the primary collisions, pN → η( ) + X, followed by 2) meson decay to scalar mediator, η( ) → π 0 S, followed by 3) mediator decay to dark matter, S → χχ. To estimate the overall meson yield, we assume 2.4 pions are produced per POT [88], while the production of η (η ) is smaller than pion production by a factor of 30 (300) [85]. We have simulated the production of η( ) mesons at the Booster with PYTHIA 8 [90]. These events are passed to a dedicated simulation which first decays the mesons to dark matter particles and then estimates the probability of passing through the detector and scattering. The MiniBooNE detector, consisting of 800 tonnnes of mineral oil, was located 491 meters downstream of the beam dump. The differential cross section for dark matter -nucleon elastic scattering in our scenario is given by where Q 2 = 2m N T , with m N the nucleon mass and T the nucleon recoil kinetic energy. The effective scalar-nucleon couplings in Eq. (12) are given by (see Eq. (A.21) in the Appendix) where in the last equality we have taken f p T u = 0.014, f n T u = 0.012 [91]. We have also assumed a dipole form factor F (Q 2 ) = (1 + Q 2 /M 2 ) −2 , with M 2 = 0.55 GeV 2 [92]. The DM search was performed in the nucleon recoil kinetic energy window of 35 MeV < T < 600 MeV. We model the detector effects with a simple step function, applying a 35% detection efficiency for events with T > 100 MeV and cutting events with smaller recoil energies. After all selections, 1465 ± 38 events were observed, while 1548 ± 198 background events were expected. With the data and background predictions in agreement, a 90% CL limit is derived by demanding the number of dark matter scattering events is less than 1.64 1548 + (198) 2 331. We have checked that this procedure reproduces well the limits derived in the vector portal DM model [87]. The parameter space in our model limited by this search is shown in the right panel of Fig. 2. We have also estimated the sensitivity of a possible future dedicated beam dump experiment using the Fermilab booster and the SBND detector [93]. The projection, also shown in Fig. 2, is made by rescaling the MiniBooNE-DM limit accounting for a factor of ∼ 20 reduction in the expected neutrino flux with a dedicated dump and the different geometric acceptance and detector mass of SBND.

C. BBN
If S is in thermal contact with the SM at the time of BBN, it contributes to the number of effective relativistic degrees of freedom at this time, which is tightly constrained by primordial element abundances. Furthermore, the decays of S to photons for m S < 2m π affects the baryon-to-photon ratio, which has the net effect of increasing the deuterium abundance D/H = (2.53 ± 0.04) × 10 −5 [94]. Recasting the bounds of [95], we find the conditions m S > 20 MeV, where the first bound is from the case where the scalar generally stays in thermal equilibrium until kinetic decoupling, while the second bound describes the limit from scalars which decay before BBN but still affect the baryon-to-photon ratio through late decays to photons.
The results of Ref. [95] do not apply if S has never reached thermal equilibrium throughout the entire history of the universe, which may occur when g u becomes very small. While a detailed analysis of the thermalization condition is rather involved (see e.g. Ref. [96]) and beyond the scope of this article, we can obtain a rough estimate from a dimensional argument: Before the QCD phase transition, the main scalar production processes are uū → S, uū → Sg and ug → uS. For small m S the S production rate is proportional to T , The scalars would reach thermal equilibrium when Γ S is greater than the Hubble rate H ∼ √ g * T 2 /M Pl . Since the number of degrees of freedom, g * , drops sharply at the QCD phase transition, we only consider T > 1 GeV. Then the thermalization condition leads to the requirement g u > ∼ 10 −9 .
For m S > 2m π and the range of couplings considered here (g u > ∼ 10 −10 ), there are no relevant bounds from BBN since S decays well before one second.
For an invisibly decaying scalar, there is in principle a bound from N eff if there are appreciable densities of S and χ during BBN. However, it is difficult to put a meaningful model-independent constraint from BBN because if S is ever in thermal equilibrium with the SM and decays to χχ, generally dark matter is overproduced in the early universe assuming standard thermal cosmology. Thus, one would have to extend the model to rectify this issue.
For instance, one could add couplings of S to electrons and/or neutrinos or introduce an extra light degree of freedom, which would provide a means to deplete the S abundance [37].
The resulting influence on BBN would then depend on the new interactions and states in the model, and so we do not show an explicit limit here.

D. Supernovae
Emission of very weakly coupled scalars could have led to increased energy loss from the core of SN 1987A, in conflict with the observed light curve. The leading contribution for the up-philic scalar S stems from the bremsstrahlung process N N → N N + S. Following Ref. [97], we treat the nucleons in a non-relativistic approximation, leading to a factorized differential cross-section where k S and E S are the 3-momentum and energy of the final-state S, respectively, and y SN N is the effective nucleon-S Yukawa coupling. Furthermore, β i,f are the non-relativistic velocities of an initial-state and a final-state nucleon, respectively. They are connected via the relationship E S = m N (β 2 i − β 2 f ). The thermally averaged energy loss rate is given by [98] where n N ≈ 1.8 × 10 38 cm −3 is the nucleon number density in the SN core, and f N is the Maxwellian nucleon distribution function. By imposing the observational bound [98] where σ abs is the thermally averaged absorption cross-section, which can be computed by again using the non-relativistic approximation.
The decay contribution is only relevant if the decay length is smaller than the size of the SN core. In our case this requires m S > 2m π , which is too massive for efficient scalar production in the SN. For completeness, we nevertheless include the decay contribution to the opacity in our numerical evaluation.
A significant fraction of the S scalars will be trapped inside the SN core if κ S is greater than the neutrino opacity κ ν [79], where κ ν ≈ 8 × 10 −17 cm 2 /g [102]. This leads to the upper edge of the exclusion region shown in Fig. 2 (top left).
The trapping bound changes for invisibly decaying scalars. Assuming that the decay S → χχ is prompt (i.e. g χ is sufficiently large), this constraint is governed by the χ + N scattering rate rather than the absorption rate of S. We will conservatively assume that a single scattering event is sufficient to trap a dark matter particle inside the SN core. Then the opacity κ χ can be computed by simply replacing σ abs in Eq. (18) with the thermal average of σ χN from Eq. (12). The resulting bound is given by the upper edge of the "SN 1987A" region shown in Fig. 2 (top right).

E. Direct detection
Conventional direct detection experiments searching for dark matter induced nuclear recoils can probe dark matter at or above the GeV scale. The effective spin-independent DM-nucleon scattering cross section is where Z (A) are the atomic (mass) number of the nuclear target, µ χN is the DM-nucleon reduced mass, and f p , f n are the effective DM -nucleon couplings, f N = g χ y SN N /m 2 S . The scalar-nucleon effective couplings were discussed above in Eq. (13) (see also Eq. (A.21) in the Appendix).
Given an assumption on the mass ratio m χ /m S and the dark coupling g χ , direct detection bounds may be shown on the m S -g u plane. In the cases where DM annihilates to SM particles, we assume g χ = 1 and choose m χ = (3/4)m S (m χ = (1/3)m S ) as a benchmark with a visibly (invisibly) decaying scalar. Then, direct detection constraints may be shown for any choices of m S and g u , as indicated on Fig. 2. There is a "thermal DM" target line where annihilation χχ → SM SM produces the observed amount of DM. Away from this region, a non-thermal mechanism is necessary to set the DM relic abundance, or if χ makes up a subdominant component of the DM (above the thermal target band), the direct detection limits which we have shown would have to be rescaled. When the relic abundance is set through secluded annihilation, choosing the scalar and DM mass immediately implies a specific DM coupling g χ , independently of g u . Therefore, we also show the direct detection limits for a secluded annihilation benchmark with m χ = 3m S . We emphasize that unlike in the visible annihilation case, the observed relic density can be obtained through secluded annihilation at every point in the m S -g u plane.
The current direct detection limits are a combination of the ν-cleus [103], CRESST [104,105], CDMSlite [106], PICO [107], and XENON1T [108] experiments. We also show the expected performance [1] of one future direct detection search aimed at observing low energy recoils, the NEWS-G experiment [109] which uses light gaseous targets. Superfluid helium detector concepts in the research and development phase may probe even lower DM masses [110][111][112].

F. Colliders
Finally, an invisibly decaying S can lead to jet plus missing energy events at the LHC through qq → Sg, S →χχ. We perform a simple parton-level recasting of the current ATLAS monojet search [113] using MadGraph 5 [114] and a modified version of a simplified model for DM coupling to a scalar mediator [115]. The resulting limit is g u < ∼ 0.1 and is independent of the scalar mass within the range we consider, as we are only considering m S much below the typical energies in the events for which ATLAS searches, on the order of hundreds of GeV of p T .

G. Summary
We summarize the limits in this section in Fig. 2 We have also looked at constraints stemming from other sources, such as neutron scattering and electroweak precision data, but found them to be not competitive with the bounds shown here.

IV. CONCLUSIONS
Light dark sectors are a particularly interesting realm of contemporary BSM phenomenology, with promising precision, beam dump, and direct detection experiments on the horizon.
The standard renormalizable portals involving a dark Higgs, dark photon or sterile neutrino have been well-examined in the literature, but much work remains to be done in identifying and investigating viable models with even broader phenomenology. In this work we have conducted a study of one such model in which a scalar mediator preferentially couples to a light SM quark, leading to a hadrophilic dark sector with a distinctive phenomenology.
We have identified the viable regions of parameter space where the dark matter abundance is set through thermal freeze-out and have performed an extensive analysis of the current experimental constraints and future probes of this scenario.
The scalar mediator we have considered is flavor-aligned but technically natural, and it represents an interesting complementary benchmark to a Higgs-like scalar. In contrast to the latter, many of the strongest bounds from rare heavy meson decays are avoided, through the absence of penguin diagrams involving top quarks with mediator emission.
Conversely, relatively large up quark couplings are allowed in our model, in contrast to the usual Higgs-like scalar mediator case where the largest reasonable Sūu coupling is of order m u /m t . This allows for significant production of scalars in light-flavor meson decays. In general, η decays at proton fixed-target experiments provide a significant source of scalars with mass below a few hundred MeV. We have examined the resulting constraints and prospects from current and future proton beam dumps, as well as precision η measurements.
At somewhat higher masses the η provides some sensitivity as well. Whether the scalar decays visibly or invisibly, relevant limits can be obtained from meson decays. Future K and η measurements at experiments such as NA62 and REDTOP, respectively, can improve these limits considerably.
In the case that the new scalar decays to SM particles, it is relatively straightforward to derive astrophysical and cosmological constraints at low m S from BBN and supernova considerations. If the scalar has a significant invisible branching fraction into dark matter, there is an inherent link between this decay channel and the dark matter relic abundance.
However, this link is highly model-dependent since the relic density may be affected by other states in the dark sector. Thus we have chosen to analyze constraints from BBN and supernovae independent of the mechanism that may set the dark matter abundance.
It would be interesting to study the relation between these in more detail in the future.
Meanwhile, at GeV-scale mediator masses direct detection starts to play an important role.
In this regard, it is interesting that future low threshold experiments such as NEWS-G offer the possibility to cover significantly wider swaths of parameter space, notably including the remaining allowed thermal dark matter target region for a visibly decaying scalar. We consider a new singlet scalar particle S. We assume the scalar interacts with quarks and gluons through the Lagrangian where the sum runs over all SM quarks ψ = u, d, s, c, b, t. This Lagrangian is defined at scales of order 100 GeV, and we consider general coupling coefficients κ ψ , κ G induced by new physics above this scale. We note that the SM Higgs couplings correspond to κ u = κ d = For the case of an up-philic scalar studied in this work, the At low scales of order 1 GeV we integrate out the heavy quarks c, b, t to obtain the effective We write these interactions in terms of the trace of the energy-momentum tensor, which is given by where we have defined In the limit of the SM Higgs couplings, we have K Θ = K u = K d = K s = 1. For the up-philic scalar considered in the main text, we have K u = (9/7)g u v/m u and K Θ = K d = K s = 0.

Matching to the chiral Lagrangian
We describe the Goldstone bosons (including the η ) with the Σ field, The pion matrix is parameterized as π = π a t a +η 0 t 0 , where t a are the usual SU (3) generators and the U (1) generator is t 0 = 1 √ 6 1. The explicit form of the pion matrix is The leading terms in the Lagrangian are capture resonant effects that become important for m S above a few hundreds of MeV. It is convenient to define the form factors π i π j | Θ µ µ |0 = Θ π (s)δ ij , π i π j | m uū u + m dd d |0 = Γ π (s)δ ij , π i π j | m uū u − m dd d |0 = Ω π (s)δ ij , π i π j | m ss s |0 = ∆ π (s)δ ij , (A.14) with which we obtain the partial width for S → ππ: where G π (s) = 2 9 K Θ Θ π + 7 9 K u + K d 2 Γ π + 7 9 K u − K d 2 Ω π + 7 9 K s ∆ π .
(A. 16) In the low momentum limit, the form factors above can be evaluated through the use of the chiral Lagrangian using Eq. (A.9). We find Θ π (s) = s + 2m 2 π , Γ π (s) = m 2 π , Ω π (s) = 0, ∆ π (s) = 0. (A. 17) We note that this agrees with the SM Higgs result in Refs. [48,51] for K Θ = K u = K d = K s = 1, but note the presence of the new form factor Ω π (s) required in the case K u = K d .
We can also compute the the partial widths for other final states. In particular, for S → K + K − and S → K 0K 0 we find where the form factors are defined analogously to those of the pions, Equations (A.14,A.16).
For ππ and KK final states the form factors at intermediate momentum scales of order 1 GeV have been extracted from data, and there are recent determinations [51] for Θ π,K , Γ π,K , ∆ π,K .
No results are available for the form factors Ω π,K , and we will therefore simply use the Chiral Lagrangian results given above. Similarly, for final states involving η, η determinations of the corresponding form factors from data are not available. At the higher scalar masses where decays to heavier mesons become important, we expect in any case that our chiral Lagrangian approach can only provide a rough estimate of the interactions of S.
In light of this, to calculate the decay width of the S at low mass we start with the results above for the ππ and KK modes, replacing the leading chiral Lagrangian form factors with those experimentally extracted from data where available [51]. In addition to these decays, we account for other modes including 4π, πη( ) and η( )η( ) with an additional channel The overall normalization of this contribution is fixed by requiring that the total S decay width to mesons matches the partonic width of S at m S = 2 GeV. Above m S = 2 GeV, we simply use the partonic width.