Dark matter relic density from conformally or disformally coupled light scalars

Thermal freeze-out is a prominent example of dark matter (DM) production mechanism in the early Universe that can yield the correct relic density of stable weakly interacting massive particles (WIMPs). At the other end of the mass scale, many popular extensions of the Standard Model predict the existence of ultra-light scalar fields. These can be coupled to matter, preferentially in a universal and shift-symmetry-preserving way. We study the impact of such conformal and disformal couplings on the relic density of WIMPs, without introducing modifications to the thermal history of the Universe. This can either result in an additional thermal contribution to the DM relic density or suppress otherwise too large abundances compared to the observed levels. In this work, we assume that the WIMPs only interact with the standard model via the light scalar portal. We use simple models of fermionic or scalar DM, although a similar discussion holds for more sophisticated scenarios, and predict that their masses should be between $\sim 100~\textrm{GeV}$ and several $\textrm{TeV}$ to comply both with the DM abundance and current bounds on the couplings of the light scalars to matter at the LHC. Future searches will tighten these bounds.

There is an overwhelming amount of evidence in favour of the gravitational existence of dark matter (DM) and its impact on baryons. Understanding the nature and origin of DM remains an open issue due to the null searches for its couplings to the Standard Model (SM) at both accelerators and astrophysically. Although purely gravitational production of DM is possible (see e.g. [1][2][3] and references therein), it is often sensitive to the assumed initial conditions in the early Universe, and it is relevant only in the absence of additional interactions that could equilibrate DM with the thermal plasma at a later stage [4].
On the other hand, the common scenario in which DM undergoes a period of thermal equilibrium with the SM species is an exceptionally attractive production mechanism. In this case, the DM relic density, Ω χ h 2 , freezes out and is essentially insensitive to the evolution of the Universe prior to the decoupling time. Instead, Ω χ h 2 is governed by an interplay between the DM mass and the coupling strength to the SM at the freeze-out tempera-ture. This scenario has given rise to a plethora of models predicting cold DM (CDM) consisting of Weakly Interacting Massive Particles (WIMPs), including e.g. supersymmetry (see [5] for a recent review) but also a number of effective models in which DM is coupled to the SM through simple portals, cf. Ref. [6].
However, the present null results of direct and indirect DM searches put ever-tightening bounds on the couplings of DM to the SM. In some of the simplest cases, the thermal value of the WIMP-like DM annihilation cross section is already excluded for masses below ∼ 100 GeV [7], while similar bounds from the Cherenkov Telescope Array [8] will reach the same level in the coming years for even larger masses, up to m χ ∼ several TeV, although this may depend on the dominant annihilation final state and the actual DM profile towards the Galactic Center (GC) [9]. The absence of any DM discovery in the near future could then point towards suppressed couplings to the SM. Unfortunately this often results in a freeze out of DM particles which happens too early in the early Universe so that then they become overproduced.
At the other end of the mass scale, ultra-light scalar fields are another type of popular species of particles beyond the SM (BSM). They appear in models related to dynamical dark energy [10] driven by quintessence [11] or emerge from modified gravity theories [12], but also e.g. in scenarios with fuzzy DM [13]. In addition, the appearance of at least one such "Stuckelberg" field could play an essential role in restoring the diffeomorphism invariance of the action in the expanding Universe [14]. It is then natural to assume that at least some of such ultralight scalars will couple to both gravity and matter and even possibly drive the accelerating expansion of the Universe [15]. Such couplings are expected to be universal to all matter species, including the SM and DM, as dictated by the weak equivalence principle, although a violation of this rule is also allowed to some extent by current observations that constrain more strongly interactions with baryons [16].
The mass of such scalars φ could be low and should be protected from possible large radiative corrections induced by φ couplings to matter thanks to a shift symmetry. A natural scenario employs the shift symmetry, φ → φ + c, where c is constant, which is preserved by derivative couplings of φ. One expects this symmetry to be only mildly broken by a non-zero mass of φ. Hence interaction terms with unsuppressed couplings, which would explicitly break the shift symmetry, are assumed to be absent.
In this study, we analyze how the thermal production of WIMP DM in the early Universe is inevitably affected by the presence of such ultra-light scalar fields coupled to both the SM and DM. To this end, we employ simple conformal and disformal shift-symmetry-preserving operators. We show that φ-mediated couplings alone can easily lead to the observed value of DM relic abundance in the regions of the parameter space of the models that are currently not excluded, and might be further testable in future searches at the Large Hadron Collider (LHC). In particular, this could provide a mechanism for models predicting too weak DM couplings to the SM and, therefore, too large abundances of WIMP DM, to be compatible with data eventually. As discussed above, such new scenarios with new scalar couplings might become of relevance in the light of current and future observational bounds.
The particular simple interaction terms that we analyze can naturally appear in the context of modified gravity, as briefly described below, but can also be considered at a more general level. We limit ourselves to the dominant operators that are normally screened from detection prospects at low energies, relevant for most terrestrial experiments and fifth force searches. On the other hand, they effectively switch on at a high energy scale corresponding to the LHC and to the early Universe when WIMP DM production occurred. This paper is organized as follows. In Sec. II, we introduce the conformal and disformal interactions of light scalars φ with matter and discuss relevant constraints, as well as basic limitations of the effective field theory approach used to study them. In Sec. III, we study WIMP DM relic density emerging from the φ-portal to the SM, while in Sec. IV we discuss its phenomenological impact on the LHC, DM, and gravitational-wave searches. We conclude in Sec. V. Selected more technical details are given in appendices. In particular, in Appendix A we discuss additional interactions that can appear in models under study, but typically play a subdominant role in determination of the DM relic density. In Appendices B and C we provide more details about solving the relevant Boltzmann equations and about the Sommerfeld enhancement factor for DM annihilations coupled to the SM via such a φ-portal, respectively. In Appendix D we discuss the scattering cross sections of disformally and conformally coupled scalars. Appendix E is devoted to discussion about a difference between gravitational and electromagnetic waves propagation in a disformally modified metric.
II. CONFORMAL AND DISFORMAL COUPLINGS OF ULTRA-LIGHT SCALARS As described above, light scalar degrees of freedom are amongst the most often discussed extensions of the SM, and they play a very special role in cosmology. A prominent example are ultra light scalar fields φ appearing in many dark energy models [10], which can be non-trivially coupled to matter [12,17]. Such interactions can naturally arise from a coupling of the φ field to the space-time metric as in scalar-tensor theories. At leading order, this affects the motion of matter fields and modifies geodesics, which can become φ dependent.

A. Light scalars and modified gravity
The most general of such coupling preserving Lorentz invariance and causality can be written [18] using the following transformation between the metric g µν in the Einstein frame, in which the gravitational part of the theory is the standard Einstein-Hilbert action of General Relativity, identified below with the Minkowski metric g µν = η µν and the one in the Jordan frame, where particles are canonically coupled tog µν , The first term in Eq. (1) is the conformal coupling while the second one proportional to D(φ, X) is referred to as disformal. The functions C and D can depend on the field φ and on the kinetic term Conformal metric transformations naturally arise in the context of scalar-tensor theories [19,20] and have been thoroughly investigated in the literature (see [17] for recent review). In particular, such theories involving kinetically dependent conformal couplings have been investigated in Ref. [21]. On the other hand, the additional disformal terms can appear in Horndeski-type [22] theories [23,24], D-brane scenarios [25,26], branon models [27][28][29] or in the context of massive gravity [30,31]. For simplicity, one often considers a purely conformal scenario, in which D = 0, or a purely disformal one with C = 1, although both types of couplings can simultaneously appear in specific models.
Light scalar fields coupled to matter can induce longrange fifth forces [11,32] which are strongly constrained for baryons [16], unless screened from local tests of gravity [15]. A fifth force in the DM sector [33][34][35] can also be constrained by observations of the anisotropies of the cosmic microwave background (CMB) radiation, large scale structures and supernovae data [36][37][38][39][40][41][42][43][44][45][46][47][48]. A large discrepancy between the fifth forces in the SM and DM sectors could also lead to observable effects related to the baryon bias [49,50]. On the other hand, non-negligible late-time interactions in the DM sector might avoid current bounds on the couplings and could be favored by cosmological data [51,52]. This could point towards interactions between the DM and a quintessential scalar field with a strength which grows in time, see e.g. [53].
At the microscopic level, interactions of φ with matter can be deduced from a series expansion of the action around the vacuum characterized by φ = 0 and the metric g µν = η µν [54,55]. The actual form of the interaction terms depend on the expansion around φ, X ≈ 0 of the functions C and D. In particular, for C(X) ≈ 1 + X/M 4 one obtains the conformal interaction Lagrangian discussed below [21], while the disformal one can be derived from e.g. exponential function D = (2/M 4 ) e βiφ 2/M 4 .
In the following, we will discuss the phenomenology of simple conformal and disformal portals between the SM and DM. While, in general, our considerations are independent of the specific theoretical motivation that lies behind these portals, it is useful to keep in mind their possible connections to other problems in contemporary physics, including the nature of the dark energy content of the Universe.

B. Models
In order to avoid model-dependent issues related to various possible UV completions of the models under study, we will follow an effective field theory (EFT) approach [56] assuming that the effects related to heavy new physics can be integrated out at the energies relevant for our discussion. We focus on the models that couple universally to all the SM species and preserve a shift symmetry of the φ field. For simplicity, we also assume that the DM particle χ is a Dirac fermion although other possibilities could also lead to similar phenomenology, and we will comment below on the case of disformally coupled scalar DM, When working within the framework of theories of modified gravity, χ DM particle moving along the geodesics corresponding to the metric given in Eq. (1) will gain φ-dependent masses. However, for small modifications introduced to the Minkowski metric, the value of m χ can be assumed constant at leading order. On the other hand, at the next-to-leading order in the m χ (φ) expansion, non-trivial couplings between χ and φ appear. We will discuss example of such interactions below, writing down effective Lagrangians that could also be studied independently of the modified-gravity context.
The DM χ particle can also have other couplings to the SM, as schematically illustrated in Fig. 1. In the following, we assume that these are subdominant, which allows us to focus on the sole impact of the φ-mediated interactions. In particular, this corresponds to the case, in which χ particles would never be in thermal equilibrium in the early Universe in the absence of the additional coupling to φ, or to the scenario in which the χ particles would equilibrate but freeze-out too early with a large relic abundance. We leave for future studies the analysis of a possible interplay between the SM-DM portals considered below and other types of couplings if both are tuned to play a comparable role. It is useful, however, to stress that even if other interactions of χ DM could cause its freeze-out with Ω χ h 2 ∼ 0.1, the impact of the coupling to φ cannot be easily dismissed from the discussion about DM relic density. Instead, as discussed below, the contribution from interactions with φ can easily be sizeable around the typical freeze-out temperatures of WIMP-like DM, T χ,fo ∼ m χ /20. derivative coupling and the disformal one where we have introduced effective conformal and disformal mass scales M SM (M DM ) that characterize the φ interactions with the SM (DM). Given our phenomenological approach, we allow the relevant scales for the SM and DM to be set independently, M SM = M DM , although our main interest is in the case when they are comparable, M SM ∼ M DM . Such universal couplings could be additionally motivated by their possible common origin. The φ field couples to the usual energy-momentum tensor of the SM and DM particles. In particular, the relevant expression for fermionic χ, which neglects terms vanishing for on-shell fermions, reads 1 where on the right hand side the indices have been symmetrized. In case of the SM, similar equation holds but with the covariant derivative, ∂ → D, which also determines the interactions between φ and the SM gauge bosons. The energy momentum tensor for a complex scalar field χ is given by On top of direct couplings to φ, additional interactions between the SM and DM particles can occur at the looplevel. However, these typically play a subdominant role when analyzing the χ relic density, similarly to possible scalar self-couplings, as discussed in Appendix A. Notably, the models described by Eqs. (4) to (7) couple pairs of φ or χ fields, therefore allowing to introduce additional discrete symmetries stabilizing the relevant particles on cosmological time scales.

C. Dominant interaction cross sections
The most relevant interaction cross sections for our discussion concern annihilation processes of a pair of Dirac fermions into a pair of light scalars,f f → φφ, where f can correspond to either the SM or DM particles. 2 In particular, assuming a negligible scalar mass, m φ m χ , the cross section for fermionic DM annihilation reads for the conformal model, where s is a square of the centerof-mass energy, and [57,58] σ ann,D,ferm. = 1 30720 π for the disformal one. Similar expressions hold for the SM fermions with M DM replaced with M SM . In particular, for both models, the annihilation cross section vanishes in the non-relativistic limit, s 4 m 2 χ , with the dominant p-wave contribution equal to where A = 32 (384) for the conformal (disformal) case, and v is the relative velocity between the DM particles. As can be seen, the thermal value of the annihilation cross section, σ ann v ∼ a few × 10 −9 GeV −2 , is obtained for a typical ratio m 6 χ /M 8 DM ∼ (10 −6 − 10 −5 ) GeV −2 assuming a characteristic WIMP velocity at freeze-out, v 2 ∼ 0.1. This ratio can naturally be obtained for m χ ∼ M DM ∼ (0.1 − a few) TeV. We will discuss the DM relic density in this scenario below in Sec. III and Appendix B. Notably, at present times χ DM annihilations into φ particles in the GC are negligible, given typical velocities of order v ∼ 10 −3 and only very mild expected impact of the Sommerfeld enhancement, cf. Appendix C.
In the discussion below, we also refer to the model with disformally coupled scalar DM χ. The relevant annihilation cross section reads [57,58] σ ann,D,c.scal. = 1 3840 π which, in the non-relativistic regime, is dominated by the s-wave contribution This cross section is larger by a factor of a few or so from σ ann v obtained for both previous models, cf. Eq. (10). Because of this, as well due to lack of velocity suppression, the thermal value of σ ann v in this model can be obtained for slightly larger values of the M DM parameter. Still, however, the correct DM relic density is obtained along the lines of constant ratio m 6 D. Current bounds and validity of EFT approach a. Collider searches The microscopic description of conformally and disformally coupled φ-mediated interactions between baryons can be used to derive constraints on such couplings based on various astrophysical probes and terrestrial experiments (see e.g. [59]). A prominent role in these efforts is played by collider searches, especially at the LHC [56,60]. The most stringent such bounds based on 13 TeV LHC data have recently been published by the ATLAS collaboration [61]. In particular, the results of the tt+E miss and disformal (cf. Eq. (5)) models. The relative weakness of the bounds for the conformal model emerge from the dominant coupling of φ to heavy top quarks and experimental challenges in performing the tt+E miss T search.
b. Validity of EFT in the SM sector These LHC constraints have been derived assuming that the validity of the EFT approach and unitarity hold in the relevant processes. In the analyses, only limited values of the centre-of-mass energy of the hard interaction have been allowed, √ŝ < g * M SM , where g * parameterizes the impact of an unknown UV completion that should not be omitted for too large momentum transfer [62,63]. In particular, when selecting events, a simplified iterative procedure to recast bounds on M SM has been applied following Ref. [64]. This is relevant for a cut-and-count analysis, while more sophisticated studies would require introducing cuts at a generator level [65]. In the following, we adopt the effective bounds from Eqs. (13) and (14) that are relevant for g * π 2 or g * π for conformal and disformal couplings, respectively [61].
We note that if g * is to be interpreted as some effective mediator coupling strength, which corresponds to the UV complete theory turning on at some high energy scale, then the aforementioned values of this parameter lie in the strongly-coupled regime. Still, however, there is a room for them not to violate the perturbativity limit, g * < 4π, which we assume to be the case in the following.
A detailed treatment of a unitarity constraint based on a partial wave expansion of the relevant amplitudes could impose additional bounds on the cut-off scales discussed above, depending on the characteristic energy of the process (see e.g. Ref [65] and references therein). In order to avoid these complications, we will limit our considerations to processes involving centre-of-mass energies √ s < O(10 TeV) characteristic for the LHC, where we assume that the EFT approach remains valid, while energies relevant for DM production will be even significantly lower.
c. Bounds on the coupling to DM Direct microscopic coupling between DM and φ remains less constrained [66]. However, in presence of sizeable such interactions in the dark sector, the fifth force between baryons could be induced at a loop level [67,68], which is tightly bounded, as discussed above. In addition, if shift-symmetry-breaking interaction terms are present, the DM-φ coupling could be further constrained by studying radiative corrections to the φ potential and its mass [69]. We note, however, that by imposing additional discrete symmetries allowed in the interaction Lagrangians Eqs. (4) and (5), the models of interest can be further protected from inducing at radiative level such shift-symmetry-breaking terms that would involve an odd number of φ fields.
In the following, we then focus on shift-symmetrypreserving conformal and disformal models with an ultralight scalar field φ and a relatively heavy DM particle, m χ ∼ TeV. Applying the EFT approach to interactions between χ and φ in these scenarios, including a particle production process in the early Universe, constrains the M DM parameter to large values, of a similar order to M SM . This suppresses any possible long-range fifth force between DM particles to be beyond current observational capabilities.
d. Validity of EFT in the DM sector In order to ensure the validity of the EFT approach in the DM-φ coupling, we perform a simple self-consistency check. In particular, as discussed in Appendix B 1, the typical freezeout temperatures of χ-φ interactions are relatively low, of order T χ,fo ∼ m χ /20, similarly to standard WIMP DM interactions. The χ DM particles then decouple being non-relativistic, and the relevant momentum transfer in the annihilation process χχ → φφ is dominated by the mass of χ. Similarly to the discussion above for the interactions between the SM and φ, we require √ s DM < g * M DM , where √ s DM is the center-of-mass energy of DM annihilations. In particular, around freezeout one obtains √ s DM 2 m χ , and the corresponding bound reads M DM c m χ , where c depends on the value of g * . For the values of g * larger than the aforementioned limits used in the ATLAS study [61], but lower than 4π, one obtains 0.1 c 1. For simplicity, in the following, we require when identifying the regime in which EFT approach remains valid, while the effects of the UV completion can be neglected. e. Other constraints relevant for ultra-high energy For large energies, additional bounds on the models under study are related to ultra high-energy (UHE) cosmicray and astrophysical probes. In particular, as shown in Appendix D, the scattering cross section of conformally or disformally coupled scalar off SM fermion can easily grow to large values for increasing energy of the incident scalar, E φ . In the rest frame of protons, the centre-ofmass energy of such collision such that √ s ∼ O(10 TeV), corresponds to E φ ∼ (10 − 100) PeV. The relevant scattering cross section can be as large as nb − µb, i.e., it can exceed the scattering cross sections of neutrinos at these energies, σ νp ∼ nb [70], by even 3 orders of magnitude. The prospects of probing such scalars produced in UHE cosmic-ray showers in the atmosphere and subsequently scattering in neutrino telescopes are, however, degraded by a small production cross section, e.g. σ(pp → φφ + X) pb for √ s ∼ 10 TeV [56]. On the other hand, both the production and scattering cross sections could be increased for even larger energy. This would, however, potentially violate EFT validity conditions and, therefore, the results of such an analysis would then be more sensitive to details of an unknown UV completion of the model.
f. Cosmological bounds on ultra-light scalars Ultralight scalars conformally or disformally coupled to the SM at temperatures around the neutrino decoupling and the onset of the Big Bang Nucleosynthesis (BBN), T BBN ∼ MeV, could contribute to the effective number of relativistic degrees of freedom and therefore change the expansion history of the Universe. This could affect successful BBN or CMB predictions. This effect is typically parameterized in terms of the number of effec-tive neutrino-like species, N eff , which is bounded from above, ∆N eff < 0.285, by the Planck data combined with the measurements of baryonic acoustic oscillations (BAO) [71].
However, as shown in Appendix B 1, for M SM 100 GeV, the decoupling of φ from the SM typically occurs at larger temperatures, T φ,fo ∼ GeV or so. As a result, the subsequent heating of the SM thermal bath around the time of the QCD phase transition generates a difference between the SM sector temperature T and the dark sector temperature of the light scalar φ. The corresponding contribution from φ to N ν at the time of BBN can then be deduced from the conservation of entropy. It is given by (see e.g. Ref. [58]) where the number of relativistic degrees of freedom at the temperature of φ decoupling T φ,fo is denoted by g * (T φ,fo ), while at the time of BBN g * (T BBN ) = 10.75. Given that g * (T φ,fo ) 60 before the QCD phase transition, one obtains ∆N eff 0.06, which is well below the aforementioned current upper limit.
Once produced in the early Universe, stable light scalar particles that decouple whilst being relativistic could contribute to a hot DM (HDM) component, on top of a similar such contribution from relic neutrinos. A detailed analysis of the allowed value of Ω φ h 2 would require marginalizing over the unknown SM neutrino masses, as well as taking into account different decoupling time and particle mass of light scalars. A simple estimate based on the upper limit on the abundance of relic neutrinos, Ω ν h 2 0.0014 [72], leads to an upper bound on the light scalar mass m φ O(eV), cf. Eq. (19) below. In the following, we will assume that φ particles are lighter than this upper limit.

III. RELIC DENSITY OF CONFORMALLY AND DISFORMALLY COUPLED WIMP DARK MATTER
As discussed above, the annihilation cross sections between the SM and light scalars, as well as between φ and DM, can be large enough so that both BSM species remain in equilibrium with the thermal bath in the early Universe. Once the temperature T drops down, the stable dark sector particles decouple from thermal plasma and freeze-out with some relic abundance, Ω χ h 2 and Ω φ h 2 , that can be obtained by solving the relevant Boltzmann equations.
Before we discuss this in more detail, it is important to note that the presence of light quintessence-like scalar fields could affect the evolution of the Universe by, e.g. introducing an early kination phase. This would then lead to an additional effect that could further modify the relic abundance of DM [73][74][75][76]. This scenario in the context of scalar-tensor theories of gravity has been first discussed in Refs [77,78] (see also Ref. [79] for further discussion) where a possible overall enhancement of Ω χ h 2 of about several orders of magnitude has been reported. Beside such a strong impact on the expansion rate at early times, the subsequent evolution of the Universe would closely resemble general relativity due to the attraction mechanism which drives the coupling between the scalar field and metric to a constant value [80][81][82]. In addition, a rapid transition to the standard cosmological scenario could further affect the DM relic abundance by initiating reannihilation processes [77]. The impact on the DM relic abundance, however, depends crucially on the coupling strength that modifies the metric. In particular, as discussed e.g. in Refs [83][84][85], stringent bounds from the Big Bang Nucleosynthesis (BBN) can introduce important constraints on such scenarios that also limit the impact on Ω χ h 2 .
In particular, the specific case of disformal coupling has recently been studied in Ref. [86], where the enhancement in the expansion rate was shown to reach up to a factor O(100). This, however, depends on the assumed initial conditions for the scalar field evolution, as well as on the value of M SM ∼ M DM , and the resulting impact can be much smaller and shifted towards large temperatures above the one corresponding to χ freeze-out. We will neglect such possible modifications to the Hubble expansion rate in the following and assume that the energy budget of the Universe around the time of χ freezeout is radiation dominated. We note, however, that, by carefully choosing the aforementioned initial conditions, Ω χ h 2 could be further changed, on top of the thermal equilibrium effects discussed below.
As discussed in more details in Appendix B, we are interested in a typical scenario where the χ DM freezeout, characterized by the temperature T χ,f.o , does not exceed the temperature of the electroweak phase transition, T EW . On the other hand, the decoupling of light scalars from the SM thermal plasma occurs at a lower temperature, T φ,f.o , although before the QCD phase transition, The last condition allows one to avoid complications due to the presence of hadronized particles among the final states of the φφ annihilations. This also helps avoiding troublesome bounds from BBN, cf. Sec. II D.

A. Boltzmann equations
The evolution of the number densities in the system composed of two dark species φ and χ/χ out of which φ is possibly much lighter and, at the leading order, the heavier species χ/χ do not couple directly to the SM, resembles the case of assisted freeze-out mechanism [87].
The relevant Boltzmann equations read then where we have simplified notations by using n χ ≡ n χ+χ . The first term on the RHS of Eq. (17) corresponds to the φφ annihilations into the SM particles, and the remaining terms are relevant for annihilations between the dark species. We leave a more detailed discussion of how to solve these equations to Appendix B, while here we briefly describe the most important features of the solutions, before discussing the results. As can be seen from Eq. (17), as long as the χ particles are in thermal equilibrium, n χ ≈ n eq χ , the φ number density undergoes a standard evolution with an additional impact from possible annihilations into χ DM that effectively add to the total thermally averaged annihilation cross section. In particular, as characteristic for any species that decouple whilst being relativistic, the yield of ultra-light φ particles in the conformal and disformal models occurs to be to roughly independent of the temperature, Y φ = n φ /s (g φ /g * s ) (45/2π 4 ), where s and g * s are the entropy density and the relevant number of relativistic degrees of freedom, while g φ = 1 is the number of degrees of freedom for φ. For φ contributing to HDM, this leads to the following relic abundance of φ where x φ = m φ /T φ,fo . Instead, for a nearly massless scalar field, the relevant present-day contribution to the energy density of the Universe is governed by the modeldependent φ potential and related slow-roll or oscillations of the field. We note that a precise determination of T φ,fo for ultralight scalars would require going to a relativistic version of Eq. (17) with the Bose-Einstein statistics taken into account. In addition, one should also include effects from a possible non-zero thermal mass of φ that is induced by scalar self-couplings (see e.g. Refs [88,89] for recent similar discussions and Appendix B for further relevant comments). However, in the analysis below, we primarily focus on the relic density of heavy WIMPs χ that freeze-out being non-relativistic. As long as n φ ≈ n eq φ close to the χ freeze-out, the relevant relic abundance is then dictated by a WIMP-like evolution employing the Maxwell-Boltzmann approximation, cf. Eq. (18), and is given by x χ,fo 20 where x χ,fo = m χ /T χ,fo , while g * (x χ,fo ) is the number of relativistic degrees of freedom at T χ,fo .

B. Results for conformal and disformal models
In Fig. 2, we show the lines of constant χ DM relic abundance, Ω χ h 2 0.12 [71], obtained by solving the Boltzmann Eqs. (17) and (18). From top to bottom, the lines correspond to disformally coupled scalar DM, conformally coupled Dirac fermion DM, and disformally coupled fermionic χ, as dictated by the relevant strength of their couplings. In each case, the correct χ DM relic density can be obtained for m 6 χ /M 8 DM ∼ (10 −6 − 10 −5 ) GeV −2 , cf. discussion in Sec. II C.
As the DM mass grows, m χ several TeV, the ratio m 6 χ /M 8 DM becomes too small unless we allow M DM m χ , which violates the approximate validity conditions of the EFT approach given in Eq. (15). This is indicated by a gray-shaded region in the plot. However, both these conditions and the χ relic abundance constraint can be simultaneously satisfied for a range of lower DM masses, m χ 3.8 TeV, 1 TeV, 300 GeV, respectively, for the three aforementioned models.
Notably, the precise value of the ultra-light scalar mass m φ , as well as the parameter M SM which governs its coupling strength to the SM, do not play important roles in obtaining the χ DM relic abundance. In particular, M SM should be sufficiently small so that φ remains in thermal equilibrium with the SM around the freeze-out of χ. This, however, is guaranteed for the range of values of this parameter above the lower limits given by Eqs. (13) and (14). As discussed in Appendix B 1, the decoupling temperature of the φ−SM interactions grows almost linearly with the increasing M SM parameter, T φ,SM ∼ M 800 GeV m χ 3.8 TeV (disformal, scalar χ). (22) In this case, as well as for a more general scenario with M DM lying below the lines with Ω χ h 2 0.12 in Fig. 2, ultra-light scalars conformally or disformally coupled to matter would have an inevitable impact on the WIMP DM relic density. In particular, the φ portal can solely provide the correct DM relic density, even if other couplings of χ to the SM are suppressed. It can also contribute to Ω χ h 2 non-negligibly in scenarios in which the other couplings are tuned to provide the DM relic abundance close to the value reported by the Planck collaboration [71]. Instead, if M DM grows to larger values, the resulting χ DM relic density obtained due to the interactions with φ remains too large, as indicated by brown and orange-shaded regions in Fig. 2. The DM abundance could then be driven to the observed value thanks to other couplings of χ, with a negligible role played by the φ portal. This also remains the case for the scenario with M DM = M SM 1.2 TeV and disformally coupled fermionic χ, where a substantial impact of the interactions with φ on the WIMP χ relic abundance would require assuming non-universal couplings with M DM about an order of magnitude lower than M SM . We note, however, that a more detailed treatment of the EFT validity limit and a possible UV completion could extend the allowed DM mass regime. A significant impact on WIMP DM relic density could then be obtained for the universal coupling to matter. This would correspond to a mass of disformally coupled fermionic DM above 2 TeV.

IV. PHENOMENOLOGY OF DERIVATIVE CONFORMAL AND DISFORMAL COUPLINGS
Given the possible important impact of conformally and disformally coupled ultra-light scalars on the DM relic density, it is relevant to discuss the phenomenological implications of such a scenario.
a. φ-SM interactions As discussed in Sec. II D, the most stringent constraints on the conformal and disformal energy scales M SM come from the null LHC searches for light scalars [56], with the dominant bounds recently published by the ATLAS collaboration [61]. These are based on 36.1 fb −1 of data collected during an initial part of the LHC Run 2. In total, one expects to collect up to ∼ 10 times more data after Run 3, and even up to a factor of ∼ 100 more data once a future High Luminosity LHC (HL-LHC) phase is taken into account [90]. Given a strong dependence of the production cross section on the conformal and disformal energy scales, σ ∼ 1/M 8 SM , this translates into an expected improvement in the relevant bounds on M SM , cf. Eqs. (13) and (14), by a factor 2.
Assuming universal couplings of φ to matter, M SM = M DM , the future LHC sensitivity to the M SM parameter could be interpreted in terms of the coupling to DM. In particular, in the full run of HL-LHC, one will be able to probe a large fraction of the mass range relevant for disformally coupled scalar DM with Ω χ h 2 0.12, cf. Eq. (22), as well as a part of the region corresponding to conformally coupled fermionic DM characterized by m χ 300 GeV, cf. Eq. (21).
In the presence of other couplings of χ to the SM, searches for φ missing energy signature could be further complemented by similar studies for χ, depending on the model. In addition, further bounds or hints of the existence of new light scalars could be deduced from the SM precision physics, e.g. the (g − 2) µ anomaly [91,92].
As discussed in Sec. II D, light scalars could also contribute to the total energy density of relativistic species in the early Universe, with a typical impact on N eff corresponding to Eq. (16) with ∆N eff 0.06. Interestingly, such a low value of ∆N eff could manifest itself in future CMB surveys, although generally at the level of only about 1σ deviation from the standard cosmological model [93]. Such an excess in N eff happens for the entire range of the M SM M DM parameter that is relevant to our discussion in Sec. III B.
b. φ-DM interactions On the other hand, searches for direct and indirect signatures of WIMP DM lie outside the high-energy and precision frontiers and, therefore, often remain unaffected by the presence of derivative conformal and disformal couplings of light scalars. In particular, as briefly discussed in Appendices A and C, both φ-loop-induced couplings and the Sommerfeld enhancement do not lead to an observable increase in the χ DM interaction rates, unless one probes regions of the parameter space of the models that lie close to, or possibly beyond, the validity regime of the EFT approach. Hence, if the χs are coupled to the SM solely via the light φ, then they remain secluded, as both their annihilation rates into visible final states and scattering cross sections off SM species are much suppressed.
On the other hand, in more general WIMP scenarios, the low-energy phenomenology can be driven by other couplings between the SM and χ DM, if they generate detectable signatures. This detection prospects, however, can also be affected by the presence of the χ DM couplings to light φs. In particular, the disformally coupled scalar DM particles have unsuppressed s-wave annihilation channels into light scalars, cf. Eq. (12). As a result, a fraction of such χs will decay invisibly, therefore reducing the expected signal rates in indirect DM searches, while the φ impact on direct detection of χ DM remains negligible. Notably, the annihilation process into φφ pairs can create a flux of boosted φs produced e.g. at the GC. This, however, for m χ ∼ TeV, remains essentially undetectable in direct searches due to tiny scattering cross sections, cf. Appendix D.
We leave a more detailed analysis of the interplay between φ-induced couplings and other portals to DM in specific WIMP models for future studies.
c. Propagation of gravitational waves An additional phenomenological aspect of the models considered in this paper appears when they are viewed as theories of modified gravity. This can be seen by writing the gravitational sector, which we assumed to be the Einstein-Hilbert action plus a scalar field φ, in terms of the standard model metric. In this frame, the theory becomes a Horndeski theory, cf. Appendix E. As a consequence, in general, gravitational waves might not propagate with the speed of light c [94]. This can lead to detectable signatures in multi-messenger astronomy that employs both gravitational-wave searches and other types of observations.
A pure conformal coupling does not generate any deviation of the gravitational wave speed from c, but a disformal coupling can have a larger impact. In particular, in the case of a pure disformal coupling to the standard model sector with D SM = 2/M 4 SM in Eq. (1), the speed of gravitational waves is given by (see Appendix E) On the other hand, the recent observations of the neutron star merger GW170817 [95] and its electromagnetic counterpart GRB170817A [96,97] constrain |c T − 1| to be smaller than 10 −15 . This can then be rewritten as a bound on the M SM parameter. Let us first assume that φ evolves slowly today witḣ φ = αH 0 M Pl , where α = O(1), H 0 ≈ 10 −42 GeV and M Pl is the reduced Planck mass. By substituting this to Eq. (23), we find a current bound M SM > 10 −8 GeV. Notably, this is a very weak constraint in comparison with the aforementioned LHC bounds and expected future reach. If, instead, the scalar field oscillates quickly around the minimum of the potential, it contributes to the DM energy density. For a harmonic potential with V = m 2 φ φ 2 /2, we have ρ kin ∼ ρ pot and hence we can estimate ρ φ ≈ ρ kin = φ2 . If the density of the φ-particles is of the order of ρ CDM , then at the present time we have φ2 0 ∼ ρ cr,0 Ω CDM ≈ 10 −47 GeV 4 . Again, by employing Eq. (23), one obtains M SM > 10 −8 GeV as it is the case for a slowly-rolling field. This is not surprising, since in both cases the time-derivative of φ/M Pl is of order H 0 . If the contribution of the scalar field to dark matter is negligible, the constraint will be even weaker.
Last, but not least, we stress that the impact of light φ remains also beyond capabilities of future multimessenger searches, since for them to become competitive with the LHC, one would require several tens of orders of magnitude improvement in the bound on the gravitational-wave speed, |c T − 1| 10 −55 .

V. CONCLUSIONS
Light new scalar φ degrees of freedom are prevalent in cosmology, as well as in various BSM models with possibly profound consequences for astrophysics. In particular, they have recently gained renewed attention in connection to the swampland conjecture [98] and the nature of the dark energy content of the Universe, while they appear also in many other contexts. If such fields exist, they could also naturally be coupled to matter.
A particularly interesting example of such a scenario employs shift-symmetry-preserving derivative couplings of φ. If such symmetry is only very softly broken for m φ = 0, then one expects that the dominant couplings of φ to other particles are induced by the aforementioned derivative operators.
In this study, we have discussed the phenomenology of such conformal and disformal couplings that are assumed to be nearly or strictly universal to all matter species. These interactions can naturally emerge from scenarios of modified gravity, in which both the SM and DM particles move along geodesics altered by φ-induced couplings. At a more general level, they can be considered as an example of higher-dimensional operators providing a portal to DM and allowing simple discrete symmetries to stabilize it.
The derivative couplings under study cause a strong energy-dependence of the interaction cross sections. As a result, the couplings effectively switch on in collider searches or at large temperatures characteristic for the early Universe, while typically remain screened from lowenergy tests. In particular, we have shown that taking into account the current bounds from the LHC and the validity regime of a simple EFT approach, the interaction rates can be strong enough to remain in thermal equilibrium around typical WIMP freeze-out. If present, they would then have an inevitable impact on χ WIMP DM relic density, with a special emphasis on scenarios with 100 GeV m χ (a few) TeV. Importantly, such sce-narios could be independently tested in collider missingenergy searches for φ, as well as they can leave a mild impact on the N eff parameter measured in future CMB surves.
We have shown that such a portal to the SM could lead to the correct DM relic density of fermionic or scalar χ on its own. Importantly, if other types of interactions between DM and the SM are also effective in the early Universe, light scalar conformal and disformal portals could still give important contributions to the total annihilation rate. The impact of φ on the present-day phenomenology of χ DM could also be non-negligible. In particular, scalar DM particles can have an efficient annihilation channel into the invisible φφ pair, which is not suppressed for small velocities. On the other hand, DM direct detection prospects would be driven by other non-φ-induced interactions.
Although light conformally or disformally coupled scalars can be well-motivated independently of models with heavy WIMP DM, both types of species can appear simultaneously in more complex scenarios. In particular, in brane-world models, disformally coupled light scalars can naturally emerge on top of heavier such species playing the role of CDM (see [58] for a recent review). On the other hand, ultra-light oscillating scalar fields with m φ ∼ 10 −22 eV could also behave like CDM themselves, as for fuzzy or scalar field DM (SFDM), see Refs. [99][100][101] for recent reviews. If such fields are coupled to matter, as discussed in our analysis, they can be active at an early epoch of the evolution of the Universe, even before rapid CDM-like oscillations begin. This could lead to a mixed scenario with both fuzzy and WIMP DM, with the latter being coupled to the SM even only through a 'fuzzy' conformal or disformal portal. We leave a detailed analysis of this scenario for future studies.
The existence of ultra-light scalar fields is particularly well-motivated in cosmology. Similarly, as far as models of DM are concerned, the dominant paradigm is to search for new heavy WIMP-like species. A possible interplay between these two extensions of the SM can have fundamental cosmological and astrophysical consequences. A more general analysis and a better connection to fundamental physics which would involve both ingredients is certainly welcome and deserves further study.  In this appendix, we briefly discuss additional interactions that can be present in the models under study, which, however, play a subdominant role in setting the relic density of WIMP-like χ DM.
a. Loop-induced couplings Besides direct interactions with the scalar field φ, the χ DM and SM particles also have loop-induced couplings to each other. However, the resulting cross sections are typically suppressed and, therefore, play a negligible role in our discussion.
In order to illustrate this, let us focus on the case of the disformal coupling that has been studied in Ref. [91]. The loop-induced effective four-point interaction between fermions in this model is given by where the explicit cut-off scale Λ has been introduced to deal with a non-renormalizable nature of the disformal coupling. One requires Λ < 4 √ π M for a perturbative treatment of the loop expansion. Assuming a universal suppression scale M SM = M DM ≡ M , one obtains from Eq. (A1) the cross section of direct χ annihilation into SM fermions, χχ → ff , of order σ dir.ann. v ∼ 10 −10 Λ 8 m 6 χ v 2 /M 16 . For Λ ∼ m χ ∼ M ∼ TeV, this yields the value of σ dir.ann. several orders of magnitude lower than the typical thermal cross section. We therefore neglect such direct annihilation processes hereafter, when discussing χ DM production in the early Universe.
The suppression is even stronger when DM annihilation in the GC are considered at present times, which is due to small DM velocities, v ∼ 10 −3 . Notably, the Sommerfeld enhancement induced by a long-range force between DM particles mediated by φ, is typically mild, S ∼ 1, cf. Appendix C. Last but not least, a loop-induced χ scattering off SM fermions, χf → χf , also remains below a measurable level, unless one focuses on regions of the parameter space with Λ > M that probe scenarios where the model-dependent impact of possible UV completions should begin to play a non-negligible role in the analysis.
A large excess of Λ over M can additionally be constrained by precision studies. A notable example is the corresponding impact on the muon anomalous magnetic moment, (g − 2) µ , which grows with increasing ratio Λ/M . On the other hand, for M 4 /Λ 3 ∼ O(1 − 10 GeV), the disformal-loop contributions might accommodate for the entire observed experimental excess [91,92]. The model would then be subject to bounds in the light of new expected data.
b. Scalar self interactions Scalar self interactions can also be induced at a loop level via virtual SM or DM particles. On top of this, if the scalar field couplings to matter are induced from an action containing φ-dependent metric, cf. Eq. (1), additional scalar self-interaction terms appear once the respective energymomentum tensor for φ is taken into account in Eqs. (4) and (5). Similarly to the couplings of φ to the SM and DM, the dominant terms are suppressed by 1/M 4 and contain terms that couple four scalar fields, e.g. (∂φ) 2 (∂φ) 2 or (∂φ) 2 φ 2 [102]. Such interactions contribute to the effective momentum transfer in the φ sector in the early Universe. On the other hand, their impact on the φ number density via possible 2 ↔ 4 processes occurs only at the next-to-leading order and is suppressed with respect to the interactions with the SM particles. We then neglect such processes below, when analyzing the relevant Boltzmann equations.
Appendix B: Details about the Boltzmann equations

Thermally averaged cross section
When solving the Boltzmann equations, an important quantity that determines both the freeze-out temperature and relic abundance is the thermally averaged cross section. In case of an approximate Maxwell-Boltzmann statistics, it can be written as [103] σ ann v T = 1 where K 1 and K 2 are modified Bessel functions of the second kind. In the relativistic limit, the proper quantum statistics should be used that can modify the collision rate.
The decoupling temperature of efficient annihilations of φ can be estimated by comparing the annihilation rate n eq φ σv ∼ T 9 fo /M 8 , which should be summed over all the relevant SM degrees of freedom, to the expansion rate of the Universe, H(T fo ) (1.67 √ g eff ) T 2 fo /M Pl . In case of relativistic scalar particles φ with a negligible mass, the typical temperatures for M SM satisfying bounds from Eqs. (13) and (14), are about GeV and slowly grow with increasing value of this parameter, T φ,fo ∼ M Pl . We have checked that this is only mildly affected when proper Bose-Einstein statistics is taken into account for φ in the annihilation cross section, cf. Ref. [88].
At the time of χ decoupling, the dominant contributions to σv for conformally and disformally coupled scalars φ annihilating into the SM species are the following: • Conformal model Due to the proportionality to the mass of the fermion in the annihilation cross section, cf. Eq. (8), the value of σv is driven by the annihilations to the heaviest kinematically available fermions, for a given center-of-mass energy √ s. The relevant expressions read where the annihilation cross section of the inverse process, ff → φφ, is given in Eq. (8). (B3) As far as χ DM annihilations into a φφ pair are concerned, the relevant annihilation cross sections are given in Eqs. (8), (9) and (11). Importantly, for m χ 100 GeV, at the time of χ freeze-out, a φ-mediated kinetic equilibrium between the SM and χ DM sectors is maintained. We, therefore, neglect the effects of a possible early kinetic decoupling when solving the Boltzmann equations.

Solving the Boltzmann equations
It is useful to rewrite the Boltzmann equations Eqs. (17) and (18) in terms of the yields of dark species, Y χ/φ = n χ/φ /s, and the x φ = m φ /T variable [87] where we have used the fact that the second term on the RHS of Eq. (B4) has to compensate for the RHS of Eq. (B5), and we have introduced In Eq. (B6), for the radiation dominated epoch, s(T ) = g * s (2π 2 /45) T 3 and H(T ) = ρ(T )/(3 M 2 P ) with ρ ρ R = g * (T ) (π 2 /30) T 4 .
For each of the dark species, the equilibrium yield is given by In particular, for a very small mass of φ we obtain i.e. the equilibrium yield remains roughly independent of T (beside a mild temperature dependence of g * s ). As a result, for a negligible mass of φ, to a good approximation Y eq φ ≈ Y φ in Eq. (B5), and the evolution of the χ DM relic density resembles the one of a standard WIMP-like scenario with the annihilation cross section dictated by the process χχ → φφ.
Appendix C: Sommerfeld enhancement for conformally and disformally coupled dark matter Light intermediate bosonic particles mediating DM self-interactions can lead to a substantial increase of the annihilation and self-scattering cross sections once they induce long-range attractive forces. This effect, known as the Sommerfeld enhancement (SE) of DM interactions [104], becomes the most prominent in the nonrelativistic limit, in which quasi-bound states can be formed.
We note that, for typical DM velocities in dwarf galaxies and the GC, v ∼ 10 −5 − 10 −3 , the potential becomes strongly coupled, V (r) ∼ m χ v 2 /2, only at very low distances, r ∼ (a few)/M , where we expect corrections going beyond the EFT approach to become important. In the following, we then assume, for simplicity, that a weakly coupled regime is valid until r r 0 = 1/M . For larger distances, the potential quickly vanishes with increasing r. It becomes effectively suppressed at a distance r φ ∼ (a few) × r 0 , which is much smaller than the typical de Broglie wavelength of a WIMP DM particle, i.e. r φ m χ v 1 for m χ M . At the quantum level, the behavior of a system of two DM particles interacting with the potential V (r) is determined by solving the relevant Schrodinger equation for the radial wavefunction. The SE factor can then be obtained by comparing the probability density of finding both particles at the same place, calculated with and without the potential, S = |ψ/ψ 0 | 2 r→0 . In Fig. 3, we show such SE factor corresponding to the = 0 term in the partial wave expansion, as a function of the DM mass. When solving the relevant Schrodinger equation, we have regularized the singular potential V ∼ r −7 by replacing it for r < r 0 with a square well with height V 0 = V (r 0 ) [106]. In a more detailed treatment, it can be shown that by adjusting V 0 as a function of the cutoff distance and requiring a smooth transition of the wavefunction between the two regimes of the potential, S can be kept independent of the choice of the cutoff [107]. Last but not least, we have normalized the SE factor so that S → 1 in the relativistic limit [108].
As can be seen in Fig. 3, for DM masses m χ M , one obtains only a very mild enhancement, S 1, as expected for the quickly vanishing potential under study. The enhancement can grow for increasing DM mass, especially close to the resonance peaks. The first such peak corresponds to m χ ∼ 3 M . This, however, already lies in the regime in which the EFT description of the microscopic interactions breaks down, and in which the potential becomes strongly coupled at larger distances, r > r 0 . A detailed treatment of such a scenario goes beyond the scope of this study. On the other hand, for the region of the parameter space of interest here, one can effectively assume that S 1.
where S SM and S DM are the action for the standard model and dark matter sectors, respectively. These sectors depend on the metricsg SM andg DM , which are related by a disformal transformation to the metric g, cf.
Eq. (1). We can write this action in terms of the metric g SM instead by employing a disformal transformation. As a result, the standard model sector is now decoupled from the scalar field φ, but the gravitational sector is no longer of the Einstein-Hilbert form. Instead, it will in general be of the Horndeski form, in which the Lagrangian reads with where X is the kinetic term for the scalar field φ. The transformation rules can be found in Ref. [23]. Under a purely disformal transformation the Einstein-Hilbert action transforms as follows, i.e. we find that in the new frame (D = 2/M 4 SM ) K(φ, X) = (1 + 2XD) 1/2 X G 3 = 0 = G 5 G 4 = (1 + 2XD) 1/2 In general, the speed of gravitational waves c T depends on the following combinations [94] F T = 2G 4 + XG 5,φ − 2XφG 5,X G T = 2G 4 − 4XG 4,X − XG 5,φ − 2HXφG 5,X with c 2 T = F T /G T . In our case we find which leads to Eq. (23) for negligible spatial field fluctuations in the vicinity of heavy objects.