Freezing-in hadrophilic dark matter at low reheating temperatures

If the reheating temperature at the end of inflation is low, of order 10 MeV, then dark matter produced through ultraviolet freeze-in has a large direct detection cross section. We study such a scenario in which dark matter is hadrophilic. This leads to dark matter-nucleon scattering cross sections of interest for near-future experiments for dark matter masses in the range of 100 keV-100 MeV. We explore how these predictions vary if reheating is non-instantaneous.


Introduction
The search for weakly interacting massive particles is mature, and direct detection experiments have placed strong bounds on this dark matter paradigm. However, flagship experiments such as Xenon-nT [1], LZ [2], and PandaX [3] have weakened sensitivity to dark matter with masses less than O(10 GeV) due to their energy thresholds. But dark matter with masses as low as O(10 keV) produced from the thermal bath of Standard Model (SM) particles can be consistent with bounds from large scale structure and Lyman-α measurements [4]. This has motivated the development of novel detector materials and technologies that can probe this interesting keV to GeV window (see e.g. [5][6][7][8][9][10][11][12][13][14][15][16][17]). Some proposals hope to probe dark matternucleon cross sections as small as 10 −45 − 10 −43 cm 2 after only ∼ 1 kg · yr exposure.
However, the number of future detection technologies may exceed the number of detectable sub-GeV dark matter candidates with concrete cosmological histories, especially those which preferentially couple to hadrons. The few sub-GeV hadrophilic dark matter candidates include asymmetric dark matter [7,18], Co-SIMPs [19], and HYPERs [20]. One challenge in constructing such models is a robust lower bound of O(10 MeV) on the mass of dark matter particles which were in chemical equilibrium with the SM bath [21]. These light particles would unavoidably contribute to the relativistic degrees of freedom at the time of big bang nucleosynthesis (BBN), altering the primordial abundances of light nuclei. Additionally, cross sections that are directly detectable for light dark matter often require light mediators. These light mediators can also contribute to dark radiation at the time of BBN.
One way around these difficulties is to consider other production mechanisms, such as freeze-in [22] wherein feeble interactions between the SM and the dark matter never result in a chemical equilibrium between the two; nevertheless, a significant dark matter density may be built up. Models utilizing this mechanism give a directdetection benchmark for sub-GeV dark matter scattering off electrons [23], and is a target for the SENSEI experiment [24].
With these future low-threshold experiments and model-building challenges in mind, we explore a model of sub-GeV hadrophilic dark matter coupled to the SM via a scalar mediator. We find that UV freeze-in [25] scenarios in which the dark matter is frozen-in instantaneously at a low reheating temperature T R through a higher dimension operator can provide interesting direct-detection benchmarks. 1 Large couplings to the Standard Model are required to reproduce the relic abundance if T R is as low as possible without disturbing the successful predictions of BBN, T R ∼ 10 MeV. 2 The predicted direct detection cross section depends on whether the maximum temperature of the SM bath during pre-heating, T max , is close to or much greater than T R . We consider in detail both instantaneous reheating and finite reheating scenarios.
The outline of this paper is as follows. In Sec. 2, we present a model of hadrophilic dark matter with a scalar mediator. In Sec. 3, we analyze the UV freeze-in of dark matter assuming instantaneous reheating and predict direct detection cross sections for future dark matter experiments for various low T R . In Sec. 4, we relax the assumption of instantaneous reheating, allowing T max > T R . Freeze-in processes with initial-state pions are especially sensitive to T max and can significantly impact the predicted cross section. We conclude in Sec. 5. In App. A, we consider variations on the model presented in the text and discuss how the choice of model impacts the predicted direct detection cross section. In App. B, we derive an expression for the contribution to the dark matter yield from T max to T R .

A model of hadrophilic dark matter
Our UV-complete model of hadrophilic dark matter contains a real scalar SM singlet φ, TeV-scale colored vectorlike fermions ψ, and a fermionic dark matter SM singlet χ. 3 The Lagrangian contains the following terms which respect a Z 2 symmetry under which ψ is odd. Because ψ is colored, it induces a coupling of φ to hadrons when integrated out, so φ is a hadrophilic mediator. The rest of ψ's SM charge assignment will impact the details of the phenomenology. For concreteness, we take ψ to transform as (3, 2, 1/6), under the SM gauge group (SU(3) c , SU(2) L , U(1) Y ), which would allow for a potential eventual embedding in a Grand Unified Theory (GUT) framework. However, we do not assume supersymmetry nor complete GUT multiplets for ψ, so substantial threshold corrections would be needed to obtain unification. For discussion how on the choice of charges for the fermions (including the possibility of adding a full GUT multiplet) impacts the results presented here, see App. A.
As long as the Z 2 breaking is small (or zero), ψ are stable over detector lengths and we may infer from the LHC searches for long-lived bottom squarks based on ionization energy loss and time of flight by the ATLAS collaboration [34].
After ψ is integrated out, φ has an effective coupling to gluons In the low-energy theory, this translates into a φ-nucleon coupling given by [35] L ⊃ −y n φnn, where m n is the mass of the neutron. Here and below, we use m p m n where m p is the proton mass. We only consider m φ much greater than the momentum transfer at direct detection experiments. Thus, the dark matter-nucleon scattering cross section relevant for direct detection is where µ χn is the dark matter-nucleon reduced mass. Since we will be focused on low T R GeV, neither the gluon coupling nor the nucleon coupling is the most relevant for freeze-in. Instead, induced couplings to pions and photons matter most. The φ coupling to pions is [35] φ also couples to other hadrons. These hadronic couplings, together with the coupling of φ to ψ (which is charged), generate a coupling of φ to photons: This coupling gets a contribution from both IR (the charged hadrons) and UV (ψ) physics. While the contribution to the coefficient of this operator from the UV is calculable, the IR contribution is subject to non-perturbative physics. For this IR contribution, we use the Naive Dimensional Analysis (NDA) [36,37] estimate ynα 4πmn . The UV contribution to this coefficient, which depends on the electroweak charge of ψ, is 15ynα 8πmn for the representation chosen above. As we will see, the induced couplings in Eqs. (2.6) and (2.7) permit dark matter to freeze-in through pion and photon annihilations.
Last, we summarize phenomenological constraints on the Yukawa couplings of Eq. (2.1) as well as m φ . For m φ lighter than the range considered below, y χ would be constrained by limits on self-interacting dark matter [38,39]. Here, however, both y ψ and y χ are only constrained by perturbativity. In our results, we set both to a maximum value of 1. The lower bound on m ψ in Eq. (2.2) implies an upper-bound on y n through Eq. (2.4), y n = 9.3 × 10 −5 y ψ 1 There are also constraints in the y n -m φ plane arising from rare meson decays. The ψ's give rise to a φtt coupling at two-loop order [7]: This coupling in turn contributes to B + → K + φ decays, K + → π + φ decays, and B −B mixing. Constraints from these processes are weaker than Eq. (2.8). The bounds from rare B + [40] and K + [41] decays have been discussed previously [7]. We have checked B −B mixing constraints [42] do not limit the allowed parameter space, and we do not discuss them further.

UV freeze-in with instantaneous reheating
Following the era of inflation, reheating could occur rapidly. If it does, the maximum temperature of the SM bath after inflation T max , could coincide with the temperature that marks the beginning of the radiation domination era T R . One cosmological history which could achieve T max T R is a multi-field inflation scenario, perhaps accomplishing reheating via a parametric resonance, see e.g. [43]. In this section, we assume that T R T max . In the next section, we discuss the implications of relaxing this assumption, allowing T max > T R .

Freeze-in calculations
In this subsection, we calculate the UV freeze-in of dark matter from SM particle annihilations for the case of instantaneous reheating. Since φ is hadrophilic, and we focus on low reheating temperatures, 5 MeV T R 20 MeV, freeze-in proceeds through photon and-perhaps less obviously-pion annihilations. For even though the abundance of pions are Boltzmann suppressed at these temperatures, they can have a significant impact. The lower bound on the reheating temperature is set by the requirement that predictions of BBN remain undisturbed [44,45]. The upper bound, as we will show below, corresponds roughly to the T R at which predicted direct detection cross sections are too small to be observable at future experiments. Note the O(TeV) ψ are not produced in the early universe due to the low T R .
We first calculate the dark matter produced from photon annihilations, γγ → χχ. The UV freeze-in proceeds through the dimension-7 operator (see discussion surrounding Eq. (2.7)), Ignoring the possible contribution from initial-state pions for the moment, the numberdensity Boltzmann equation governing dark matter freeze-in is [46] n χ + 3Hn χ ≈ σ χχ→γγ v n eq χ 2 , σ χχ→γγ v n eq where |M| 2 is the fully averaged matrix element squared for the process. The 1/2! appears due to the 2 identical photons. In the second line, g χ and g γ are the degrees of freedom of dark matter and photons. The phase space has been simplified to a single integral since |M| 2 is only a function of s (see e.g. [25,47]). This assumes Maxwell-Boltzmann distributions for thermal γ's and χ's and that χ's are always out of equilibrium. Proceeding in a similar vein, we can calculate the UV freeze-in contribution coming from pion annihilations π + π − →χχ due to the couplings in Eq. (2.6). Adding both contributions and integrating over temperatures, we obtain the yield of dark matter, Y DM = n DM /s, the ratio of the dark matter number density to the SM bath entropy. The total yield obtained from pions and photons is where κ j = 1 for j = π ± , 17α 4π for j = γ, (3.4) and Above, M Pl = 2.4 × 10 18 GeV is the reduced Planck mass, g (s,) * are the (entropic) relativistic degrees of freedom, and Γ φ is the total width of φ which in our case almost always decays to χχ. The presence of this width as well as s/m 2 φ in the denominator of Eq.  where Ω DM is the energy fraction of dark matter, ρ crit is the critical energy density, and s 0 is the entropy of the photon bath today [48]. For the low reheating temperatures we consider, one might suspect that photon annihilations are the dominant process for dark matter freeze-in. For T R 13 MeV, this is the case, as shown in Fig. 1. This figure displays the yield of dark matter frozen-in via pions divided by the corresponding yield frozen-in via photons as a function of the reheating temperature, i.e., the ratio of the two terms in the sum of Eq. (3.3). The different colored curves correspond to different possible SM charge assignments for ψ. The quark doublet-like representation we have emphasized is displayed in green. The band around the curves corresponds to varying the IR contribution to the photon coupling around the NDA estimate of Eq. (2.7) by a factor of 2. As expected, the width of the band is smaller for those representations that receive larger contributions to φF F from integrating out ψ in the UV. For now we continue to concentrate on ψ ∼ (3, 2, 1/6), and we postpone discussion of other representations to App. A.
Returning to Fig. 1, we note that for temperatures T R 13 MeV, pion annihilations are the dominant freeze-in channel. 5 This is despite the exponential penalty required to find two pions in the relatively cold thermal bath. The presence of this additional channel reduces the predicted direct-detection cross sections since it requires smaller couplings than would naively be predicted if only photon annihilations were taken into account. It also foreshadows the importance of the assumption T max = T R that we have made thus far. Allowing for T max > T R can substantially increase the yield from pions, as we will see in the next section.
One final comment: electrons will also be in the SM bath at low T R and in principle could contribute to the freeze-in of dark matter. However, relative to the yield from photon annihilations in Eq. (3.3), we expect this yield to be suppressed by at least α 2 m 2 e /T 2 R and is therefore negligible.

Results
For fixed m χ and T R , explaining the dark matter abundance in Eq. (3.6) predicts a value of y n /m 2 φ . We plot this requisite value as a function of m φ in Fig. 2 for a benchmark case of (m χ = 0.5 MeV, y ψ = 1, y χ = 1) and various reheating temperatures. The shaded bands surrounding these curves (which are noticeable only for lower reheating temperatures) correspond to varying the NDA estimate for the photon coupling in Eq. (2.7) by a factor of 2. Uncertainty in the photon coupling due to NDA estimates changes the predicted σ χn by roughly 15% for lower reheating temperatures, and this effect gets smaller as the contribution from the pions outstrips that of the photons at higher temperatures. We will not show this uncertainty in the rest of our results.
The constraints on y n and m φ discussed in Sec. 2 are shown in shades of gray. The light gray shaded region is ruled out by current direct searches for ψ from the LHC; see Eq. (2.8). The dashed gray line assumes a possible future bound from the LHC of m ψ 2 TeV. We also show the approximate lower bound on T R coming from BBN in light gray. The darker gray regions are constrained by searches for rare B and K decays. These constraints are subdominant except for a narrow range of φ masses near m φ ∼ 200 MeV.
A noteworthy feature of this plot is the resonance around m φ ∼ 2m π + . This resonance causes pion annihilations to become more productive and diminishes the requisite nucleon coupling and cross section. Because one point of interest is to find benchmarks for future direct detection efforts, to achieve larger cross sections and simplify our analysis, we require m φ ≥ 1 GeV in the rest of this section. For these heavier φ's, the requisite value of y n /m 2 φ is roughly independent of m φ (as evident in Fig. 2). In our numerical computations, we do not assume UV freeze-in and  calculate the freeze-in with the full φ propagator, allowing for its on-shell production and subsequent decay to pairs of dark matter. However, once we impose this lower bound on m φ , the effective operator picture of Eq. (3.1) is valid [49] and allows our results to conform to some of the usual intuition from simpler UV freeze-in scenarios.
The right vertical axis of Fig. 2 shows that low reheating temperatures may predict detectable cross sections. Going to large T R rapidly decreases the cross section, even more quickly than if only photon annihilations were considered. Despite their exponentially small abundance in the thermal bath, pions quickly overtake photons to become the dominant driver of freeze-in for T R 14 MeV (see Fig. 1). Nevertheless, it is informative to first consider the limiting case in which T R 12 MeV so that pions contribute less than ∼ 1% to the dark matter relic abundance and may be safely ignored. Then the total dark matter yield comes from photon annihilations alone and for m χ T R , the integrals in Eq. (3.3) may be done analytically, yielding:  Thus, lower reheating temperatures predict potentially detectable cross sections. The full predictions, including the pion contributions, for the direct detection cross section as a function of m χ are shown in Fig. 3. The colored contours correspond to the same reheating temperatures as in Fig. 2, the gray regions correspond to the same bounds, and y ψ = y χ = 1 again. The LHC constraint again comes from direct searches for ψ (see Eq. (2.8)), and the dashed gray line assumes instead a possible future bound of m ψ 2 TeV from the LHC.
Note that the approximation of Eq. (3.8) holds as expected when T R 12 MeV and m χ T R . As m χ increases, m χ T R , the photons able to produce the χs start becoming Boltzmann suppressed, requiring larger couplings and resulting in a larger predicted cross section. Also note the increased spacing between contours for T R 14 MeV, when pion annihilations start to dominate.
Much of the viable parameter space corresponding to 6 MeV T R 16 MeV will be testable by a host of proposed direct detection experiments whose projected sensitivities are shown with dashed gray lines [5, 10, 12-14, 16, 17]. We emphasize the ability of this cosmological history to predict cross sections which will be detectable in the future for dark matter masses below O(10 MeV). In fact, dark matter may be as light as 200 keV and still be detectable by a 1 kg · yr exposure of CsI [13].
In Fig. 3, we have chosen m φ ≥ 1 GeV for simplicity. As shown in Fig. 2, the LHC bound stops dominating for m φ 1 GeV. For smaller φ masses, the LHC constraint below the BBN constraint in Fig. 3 would be absent. This would allow dark matter masses as light as ∼ 80 keV to be probed by future detectors if T RH ∼ 6 MeV. We do not present this possibility in detail, however. For such light φ masses, the contours for T R 8 MeV are sensitive to m φ due to the pion resonance in Fig. 2.

UV freeze-in with finite reheating
In the previous section, the total dark matter yield was obtained assuming that reheating happened instantaneously at the temperature T R . However, reheating from perturbative inflaton decay may not be immediate; 6 then the yield in Eq. (3.3) can significantly underestimate the actual dark matter yield [52,53]. Indeed T max may be much larger than T R ; assuming a matter dominated era of pre-heating T max 0.5 (m Φ /Γ Φ ) 1/4 T R (see e.g. [53]). Since 10-MeV-scale reheating corresponds to Γ Φ ∼ O (10 −23 ) GeV, a simple single-field inflation model with instantaneous reheating from perturbative decay achieving T max ∼ T R would imply m Φ ∼ Γ Φ is ultra-light. Having T max > T R may allow additional inflationary scenarios and also open new possibilities for baryogenesis. In this section, we redo the dark matter yield calculation in the case T max > T R and detail the predictions for direct detection.

Freeze-in calculations
While a significant amount of dark matter can be produced at temperatures above T R , it is also diluted by the subsequent entropy generated from inflaton decays. This interplay between enhanced production and entropy dilution means that the yield of dark matter generated between T max to T R is sensitive to the temperature dependence of the freeze-in process(es) [53]. We find the contribution to the dark matter yield  from T max to T R is where the I(m j ) is defined in Eq. (3.5). For the pion contribution, T max is capped at the temperature of the QCD phase transition, 156.5 MeV [54]. Above this temperature, gluon annihilations to dark matter occur, but we find these contribute negligibly to the final dark matter relic abundance. The pre-factor of 0.4 is a correction that comes from numerically solving the Boltzmann equations [53]. Combining this yield with the one from temperatures below T R , found in Eq. (3.3), gives the total yield of dark matter in inflationary scenarios with arbitrary T max . For the small values of T R we consider, a larger T max implies exponentially more pions in the SM bath. The resulting temperature dependence of the thermallyaveraged pion annihilation cross section can more than compensate for the entropy dilution between T max and T R . Thus, the dark matter yield can be much larger than what we found assuming T max T R .
To illustrate this point, in Fig. 4 we show the ratio of the total dark matter The solid lines in Fig. 4 prove that the total yield can be many orders of magnitude larger than expected from just the instantaneous reheating calculation, especially if T max is appreciably larger than T R . As expected, this is due to the greatly enhanced contribution from pion annihilations, as is evident from the dashed contours. The dotted lines show that the contribution to the dark matter yield from photon annihilations from T max to T R is O(1) times its contribution from T R and below. Interestingly, the difference between the solid lines with T max /T R = 5 and ∞ is much less than between T max /T R = 5 and 2. This, together with the previous observations, emphasizes that for large T max /T R , the dark matter relic abundance is set by pion annihilations above T R . Furthermore, the contributions from the photon annihilations are always subdominant and can be ignored when T max 2T R .

Results
As in Sec. 3.2, we plot the value of y n /m 2 φ which gives the correct relic abundance as a function of m φ . In Fig. 5, this is done for a benchmark case of (m χ = 10 MeV, y ψ = 1, y χ = 1) and various reheating temperatures. The bands for each T R value correspond to 5 ≤ T max /T R ≤ ∞, with the top of each band 5 and the bottom of each band ∞. Also shown in shades of gray are the same constraints from Fig. 2 and discussed in Sec. 2.
It is clear that the sensitivity to the pion resonance increases for larger T max /T R . We limit our parameter space scan by imposing m φ ≥ (2, 3) GeV for T max /T R = (5, ∞), respectively, for the rest of this section. This both achieves larger cross sections and simplifies the analysis. As is evident from the figure, the requisite value of y n /m 2 φ is then roughly independent of m φ . The right vertical axis of Fig. 5 shows that direct detection cross sections, while still detectable for lower T R , are greatly suppressed as T max /T R increases, a result of the enhanced yield from pions.
The full predictions, including the dominant contributions from T max to T R , for the direct detection cross section as a function of m χ are shown in Fig. 6 for both T max /T R = 5 (in the top panel) and T max /T R = ∞ (in the bottom panel). The colored contours, constraints, and future-experiment sensitivities correspond to those in Fig. 3. Again, we have set y ψ = y χ = 1. Notably, contours of fixed T R now occur at much lower values of y n /m 2 φ , and so the LHC constraint constrains a smaller area of the space allowed by BBN.
As expected, the effect of the enhanced freeze-in from T max to T R due to pion annihilations is to suppress the predicted direct detection cross section. The result is that future experiments will only be sensitive to a small range of parameter space for these inflationary histories. Still, there are benchmarks that predict detectable dark matter as light as ∼ 1 MeV if T max /T R = 5 or ∼ 2 MeV if T max /T R = ∞. Note also, for any inflationary scenario which results in a given T R , the predicted cross section will be greater than or equal to the one predicted by T max /T R = ∞.

Discussion
We have explored a model of hadrophilic, sub-GeV dark matter with cross sections detectable at proposed direct detection experiments. After taking into account all collider, BBN, and rare-meson-decay bounds, we find that UV freeze-in at O(10 MeV) reheating temperatures predicts sizable interactions for dark matter as light as 200 keV. This is significant because only a few dark matter models exist at such light masses which 1) may be detected at these proposed experiments and 2) are not already ruled out by cosmological or other constraints. Some proposals for detecting light dark matter have even resulted in experimental collaborations, e.g. using superfluid helium as a target [55,56]. The existence of consistent benchmarks are important to support the science cases of these experiments.
We have determined that even at such low reheating temperatures, contributions from pion annihilations can dominate the relic abundance. Thus, we have also shown that if the SM bath reaches a higher temperature prior to reheating, T max > T R , the corresponding cross sections are reduced. Even still, these more generic inflationary scenarios can predict benchmarks of relevance for direct detection proposals.
Given the generic challenges associated with constructing light dark matter models which are within the reach of the next generation of experiments, our lowreheating UV freeze-in scenario through a scalar mediator represents a significant benchmark for these proposals. For instance, two simple anomaly-free vector-mediator models, a kinetically mixed dark photon and gauged B − L, result in cross section many orders of magnitude below future sensitivities. This is due to the relative ease with which electron-positron pairs annihilate to produce dark matter through a dimension-six operator, in contrast with our scalar-mediator model.
In this work we found the enhanced UV freeze-in production from non-instantaneous inflaton decays to have a significant effect on the expected direct detection cross section. It would be interesting to consider dark matter production from UV freeze-in (and associated phenomenology) in more exotic models of reheating. 7 Here, for simplicity, we have assumed that Z 2 breaking is absent. However, this need not be the case, and if present, would potentially allow Yukawa couplings between Standard Model fermions and the ψ fields. Such couplings could allow ψ to impact precision electroweak observables, and could also change the details of the bounds that come, e.g. from meson decays, or direct LHC searches. It would be interesting to consider the phenomenology of Z 2 -breaking scenarios in more detail.

A Alternative UV completions of hadrophilic dark matter
In this appendix, we consider alternative UV completions of the hadrophilic dark matter model. In particular, we consider having one or more vectorlike pairs of quarks and leptons, ψ i + ψ i , in various representations under the SM gauge group, and we will study how the phenomenology depends upon the choice of representation. The generalized version of the Lagrangian in Eq. (2.1) includes The induced coupling of φ to gluons, after integrating out the colored ψ i + ψ i , is now given by: The coupling to gluons can then be mapped onto the low-energy nucleon coupling y n via 8 : Electromagnetically charged ψ i also induce a one-loop coupling of φ to photons: where q(ψ i ) is the electromagnetic charge matrix of ψ i in the weak-isospin space, and d c (ψ i ) is the dimension of ψ i in SU(3) c color space. Due to the coupling to gluons in Eq. (A.2), φ also has induced couplings to other charged hadrons, which together with the coupling to protons in Eq. (A.3), generate a contribution to the φF µν F µν operator in the IR-this is in addition to the contribution from integrating out ψ i 's that carry an electromagnetic charge. Therefore, the coupling of φ to photons can be parameterized in terms of the nucleon coupling, as: where c IR γ = 1 (NDA estimate), (A.6) 8 More generally, if 1/Λ g is the coefficient of the φG a µν G µν,a operator, then the low-energy nucleon coupling is given by y n = 8πmn b0αsΛg . Here b 0 = 11 − 2 3 n L is the coefficient of the QCD β-function at one-loop with n L being the number of quarks lighter than Λ QCD and m φ [35]. , and y ψ i /m ψ i of various ψ i . The total dark matter yield obtained from pions and photons from temperatures below T R (from T max to T R ) can be computed using Eq. (3.3) (Eq. (4.1)) with κ π ± = 1 for the pions and κ γ = α 2π c IR γ + c UV γ , generalizing the contribution from photons in Eq. (3.4). As discussed in the main text, the photon contributions are negligible compared to the pion contributions for large T R (T max ). In the case where T max T R (i.e. instantaneous reheating), for T R 12 MeV where the pion contributions can be ignored σ χn is 17 UV completion with a single pair of ψ + ψ: If we restrict ourselves to the case where there is only one vectorlike quark pair ψ + ψ that couples to φ as in Eq. (2.1), then the UV contribution to φF µν F µν operator is given by independent of y ψ /m ψ . In the main text, ψ is taken to transform as SM-like weakisodoublet Q ∼ (3, 2, 1 6 ) under the SM gauge group for which c UV γ = 15/2. For other SM-like representations for ψ, namely up-type weak-isosinglet U ∼ (3, 1, 2 3 ) and down-type weak-isosinglet D ∼ (3, 1, − 1 3 ), c UV γ = 12 and 3, respectively. Since, by assumption, the vectorlike quarks in these UV completions do not mix with the SM quarks, they can be long-lived and stable over detector lengths. Therefore, m ψ for various representations of ψ can be constrained from the searches for pair-production of quasi-stable vectorlike quarks at particle colliders. In particular, we can infer bounds on up-type (down-type) vectorlike quarks from the LHC searches for long-lived top (bottom) squarks, by the ATLAS collaboration [34], to be If reheating is assumed to be instantaneous (i.e. T max T R ), it is evident from Fig. 1 that the specific representation of ψ somewhat matters only when T R is less than around 12 MeV, as the pion annihilations start dominating for larger T R . In comparison with the results in Fig. 3, the colored contours for T R 12 MeV move slightly (down) up if ψ transforms as (U ) D instead of Q. Specifically, for m χ T R , the cross sections σ χn that yield the correct relic abundance for UV completions with U + U and D + D are (17/26) 2 and (17/8) 2 times the corresponding cross sections for Q+Q in Fig. 3. The bounds on σ χn for up-type and down-type weak-isosinglets, based on Eq. (A.10), become stronger by a factor of 4 compared to the exclusion region labeled as LHC in Fig. 3 for the SM-like weak-isodoublet quark representation.
Finally, we also comment on what happens if ψ has more exotic quantum numbers: ψ 0 ∼ (3, 1, 0). Since ψ 0 is an electroweak-singlet, c UV γ = 0 from Eq. (A.7), and therefore the φF µν F µν operator is solely generated from IR physics (i.e. charged hadrons). As a result, the dark matter yield from photon annihilations becomes more sensitive to the uncertainty in the NDA estimate of c IR γ = 1. In this case, hadronization will result in fractionally charged hadrons which are not presently the target of a dedicated search at ATLAS/CMS. And while limits relying on dE/dx do not immediately apply, searches looking for slow moving particles using time of flight measurements should. Based on the ATLAS searches for long-lived R-hadrons using the time of flight measurements only [34], we estimate: m ψ 1.5 to 1.7 TeV (95% CL; electroweak-singlet ψ), where the range of values comes from assuming a signal selection efficiency that is a factor of 2 higher or lower than that of the long-lived top-squark searches. Here, we obtained the observed signal upper limits using the Zstats package [59] based on Bayesian-motivated statistical measures in [60]. A dedicated search might extend these limits slightly. Note that σ χn with ψ ∼ ψ 0 , for m χ T R and T R 12 MeV, is enhanced by a factor of (17/2) 2 compared to the corresponding contours in Fig. 3. The bound on σ χn using the y n bound of Eq. (A.10) and m ψ 1.7 TeV is roughly a factor of 4 stronger than the corresponding exclusion labeled as LHC in Fig. 3.

UV completion with a full GUT multiplet of vectorlike fermions:
Finally, we consider UV completions with vectorlike fermions in 5 + 5 and 10 + 10 representations of SU(5), where 5 = D + L and 10 = Q + U + E contain all of SM quark and lepton representations. Here, L + L and E + E are SU(2) L doublet and (charged) singlet vectorlike leptons, respectively. This possibility might be motivated for possible embedding in a GUT framework, and could preserve the apparent gauge coupling unification if we assume supersymmetry. We, however, do not assume an unbroken SU(5) gauge symmetry in the UV.
In the special case where y ψ i /m ψ i are the same for all vectorlike fermions, we obtain the UV contribution to φF µν F µν operator at low energy to be: For vectorlike fermions in 5 + 5 and/or 10 + 10 of SU(5) with the same y ψ i /m ψ i , c UV γ = 12. In a realistic GUT scenario, however, y ψ i /m ψ i for various vectorlike fermions in a GUT multiplet at low energies, after the renormalization group evolution from the GUT scale, would not be same. After specifying a GUT, it is straightforward to obtain c UV γ in terms of I c , d c , d L , Tr[q 2 ], and y ψ i /m ψ i of various ψ i . Fig. 7 shows the predicted direct detection cross section σ χn as a function of the dark matter mass m χ for heavy vectorlike fermions in complete SU(5) representations: 10 + 10 in the top panel and 5 + 5 in the bottom panel. For simplicity and illustration purposes, absent a concrete GUT model, we chose y ψ i /m ψ i (≡ y SU(5) /m SU(5) ) to be the equal at the IR scale for all vectorlike content in a chosen SU(5) pair and set y ψ i = y χ = 1. The colored contours correspond to various reheating temperatures from 6 to 18 MeV. The constraints and future experimental sensitivities labeled in the plot correspond to the ones shown in Fig. 3. The constraint labeled as LHC, in particular, is obtained from an analog of Equation (A.10), with d L (ψ)/2 replaced by i I c (ψ i )d L (ψ i ), with y SU(5) 1 constrained only by perturbativity and m SU (5) 1.5 (2) TeV for shaded region (long dashed line). Under these assumptions, we can note the slight increase (decrease) in the parameter space for 10 + 10 (5 + 5) of SU (5)

B Finite Reheating Yield Calculation
This Appendix derives the yield of dark matter due to a period of finite reheating found in Eq. (4.1). We make use of the notation and derivations in [53]. There are two major changes in this non-instantaneous reheating derivation relative to the instantaneous one: 1) the usual (approximate) ∂T /∂t = −HT no longer holds and 2) Y χ = n χ /s is no longer a useful parameterization since bath entropy is not conserved as the inflaton decays. Updating the first relation, we find that while the inflaton is decaying: Γ Φ is the inflaton decay rate and the "end" subscript corresponds to the end of inflation. We have in mind a single-field inflationary scenario and assume that at the end of inflation, the inflaton undergoes coherent oscillations about its potential minimum.
As for the second change, it is useful to note that: where the left hand side is evaluated at the reheating temperature. Assuming that the inflaton dominates the energy density of the Universe until the end of reheating allows us to rewrite a a end The O(1) number depends on the particular inflationary model, but drops out from the final yield expression. We further assume that A v 1 to find that