A Global View of the Off-Shell Higgs Portal

We study for the first time the collider reach on the derivative Higgs portal, the leading effective interaction that couples a pseudo Nambu-Goldstone boson (pNGB) scalar Dark Matter to the Standard Model. We focus on Dark Matter pair production through an off-shell Higgs boson, which is analyzed in the vector boson fusion channel. A variety of future high-energy lepton colliders as well as hadron colliders are considered, including CLIC, a muon collider, the High-Luminosity and High-Energy versions of the LHC, and FCC-hh. Implications on the parameter space of pNGB Dark Matter are discussed. In addition, we give improved and extended results for the collider reach on the marginal Higgs portal, under the assumption that the new scalars escape the detector, as motivated by a variety of beyond the Standard Model scenarios.


Introduction
The hypothesis that the Dark Matter (DM) is composed of weakly interacting massive particles (WIMPs) whose abundance is determined by thermal freeze-out has been a leading paradigm for decades. The direct searches for DM scattering on nuclei, however, have kept reporting null results, lowering the cross section limits at an impressive pace and ruling out many WIMP models along the way. Currently, the strongest sensitivity has been achieved by the XENON1T experiment [1].
A compelling exception to this tension is obtained if the WIMP arises as a pseudo Nambu-Goldstone boson (pNGB). This possibility is especially motivated by theories where the Higgs and the DM arise together as composite pNGBs, thus addressing simultaneously the naturalness and DM puzzles [2]. Models of this type have received increasing attention lately [3][4][5][6][7][8][9][10][11][12][13][14][15][16]. The key point is that the scalar DM, assumed to be a singlet under the Standard Model (SM) gauge interactions, couples to the SM through higher-dimensional, derivative interactions with the Higgs field, that arise from the kinetic term of the nonlinear sigma model. Up to dimension six the relevant Lagrangian contains the derivative Higgs portal, where φ is the scalar DM candidate and H the Higgs doublet. The scale f is the common decay constant of the pNGBs, while c d is an O(1) coefficient. For definiteness we have taken φ to be a real scalar stabilized by a Z 2 symmetry. The important observation is that the operator proportional to c d mediates s-wave annihilation to SM particles, but is very suppressed in the scattering of DM on nuclei, which is characterized by low momentum transfer |q| 100 MeV. Hence, thermal freezeout can yield the observed DM density for m φ ∼ O(100) GeV and f ∼ O(1) TeV, while the direct detection signal is beyond the foreseeable experimental reach. In complete models additional terms can be present beyond Eq. (1), which either respect or break explicitly the DM shift symmetry (some explicit breaking is of course necessary, to generate a nonvanishing m φ ), but the derivative Higgs portal is the irreducible ingredient underlying this type of DM candidate. A realization of pNGB DM can also be obtained [17] by adding to the SM a complex scalar field with potential invariant under a U (1), which is softly broken only by a mass term (notice that in this setup the Higgs field is not a pNGB). The angular mode of the scalar is stabilized by a remnant Z 2 and is identified with the pNGB DM. As we briefly discuss in Sec. 2, if the radial mode is sufficiently heavy it can be integrated out, yielding Eq. (1) as low-energy theory. See Refs. [18][19][20][21][22][23][24][25][26] for recent studies of models built following this approach.
Experimental probes of the derivative Higgs portal include indirect DM searches in cosmic rays, and collider experiments. In this work we explore the latter for the first time, by considering the production of two invisible φ particles through an s-channel Higgs boson. 1 While for m φ < m h /2 the sensitivity on f /c 1/2 d is immediately obtained from the available studies of invisible Higgs decays, for m φ > m h /2 the Higgs is off-shell and the momentum dependence of the coupling in Eq. (1) has an important impact on the signal kinematic distributions, mandating a dedicated study that is the main subject of this paper. We focus on Higgs production via vector boson fusion (VBF), see Fig. 1, because it is expected to provide the best sensitivity both at high-energy lepton colliders and at hadron colliders. We will expand on this momentarily.
Throughout our analysis we also consider the marginal Higgs portal, which provides an important term of comparison. 2 We assume that φ escapes the detector and thus manifests as missing momentum. While the hypothesis that φ is a thermal relic that interacts with the SM through Eq. (2) is mostly ruled out by direct detection (exceptions being the Higgs resonance region, or DM heavier than a few TeV), non-standard but motivated cosmological histories can open up large regions of parameter space [28]. Furthermore, new scalars with sizable couplings to the Higgs are broadly motivated by open problems of the SM other than DM, such as baryogenesis and naturalness. These scalars may be stable on collider timescales or decay invisibly, giving rise to the signature studied here. This scenario was discussed in the context of electroweak baryogenesis in Ref. [29]. Furthermore, natural models where the top partners are scalars require that they couple to the Higgs with strength fixed by y 2 t . If the scalar top partners are neutral under all SM gauge symmetries, as in the recently-proposed Neutral Naturalness theories of Refs. [30,31], then probing the interactions L − y 2 t |H| 2 (|ũ c 1 | 2 + |ũ c 2 | 2 ) through off-shell Higgs could be key to discovering this type of solution to the little hierarchy problem.
Previous studies of the off-shell marginal Higgs portal to invisible scalars are available in the literature. The sensitivity on λ at lepton colliders was carefully analyzed in Ref. [32] for √ s ≤ 1 TeV, where Zh associated production mostly dominates (see also Refs. [33,34] for related studies), and in Ref. [35] for √ s = 1, 5 TeV, considering ZZ fusion production. 3 The sensitivity at the LHC and FCC-hh was thoroughly examined in Ref. [36], including the VBF, monojet and tth production modes.
Here we focus on VBF production, which at lepton colliders provides the leading sensitivity for √ s 1 TeV, and at hadron colliders was found to be superior to monojet and tth in the previous study of Ref. [36]. We provide the first results for the derivative Higgs portal, in the form of sensitivity projections for a wide set of colliders that include both high-energy lepton machines and hadron machines. The collider parameters we have assumed are summarized in Table 1 The excluded region is the one enclosed by each line. All limits are derived from VBF. Right: marginal Higgs portal. The excluded region is the one above each line. Solid lines correspond to VBF processes, dashed lines to Zh associated production. The dotted line labeled "Naturalness" corresponds to λ = √ 4N c y 2 t , namely two degenerate scalar top partners.
LHC. Our VBF analysis is presented in Sec. 3, but we anticipate its main results in Fig. 2, which provides an overview of the sensitivity on both types of portals. The reach along benchmark slices of the parameter space is summarized in Table 2. In the presence of N real scalars, degenerate and with identical couplings, their effect on the signal studied in this paper is equivalent to that of a single field with {c d , λ} → √ N {c d , λ}. For example, for complex pNGB DM χ with derivative portal ∂ µ |χ| 2 ∂ µ |H| 2 /f 2 ∈ L we should take an effective interaction strength c d = √ 2 , whereas for degenerate scalar top partners (two complex scalars in the SU (N c ) fundamental, with y 2 t coupling to the Higgs) we have λ = √ 4N c y 2 t . Before concluding this Introduction, we briefly summarize the reach for m φ < m h /2, in which case on-shell invisible Higgs decays are the main search avenue. These have been analyzed in great detail in the literature, and the resulting bounds on the Higgs branching ratio to new invisible particles 4 are reported in the upper row of Table 3. At the LHC the strongest sensitivity is obtained in VBF production [37]. A future Higgs factory will increase the reach significantly, by exploiting Zh production, and an even stronger bound can be achieved at FCC-hh. The analysis of on-shell decays is insensitive to the specific portal considered, and the limits on BR(h → inv) are straightforwardly translated into constraints on f /c 1/2 d and λ by using the expression of the width, The results are given in the middle and bottom rows of Table 3. The remainder of this paper is organized as follows. In Sec. 2 we shortly discuss the pNGB DM framework, using a general effective field theory (EFT) for the SM plus the scalar DM. Section 3 presents our analysis of the off-shell derivative and marginal Higgs portals, exploiting the VBF channel at future colliders. We discuss our results and the outlook in Sec. 4. Appendix A provides details on the analysis, including a summary of all the selection cuts we imposed and an explicit comparison with the results of Ref. [36]. Finally, Appendix B gives general cross sections for pNGB DM scattering on nuclei and annihilation.

Pseudo Nambu-Goldstone boson Dark Matter
We focus on theories where both the Higgs doublet and a scalar DM candidate, singlet under the SM gauge symmetries, arise as light pNGBs from a strongly coupled sector. In this section we take the DM to be a complex scalar χ ; for real DM φ , one simply replaces χ → φ/ √ 2 . We consider the effective Lagrangian where the DM-SM interactions are described in general, up to dimension six, by Here χ * ↔ ∂ µ χ ≡ χ * ∂ µ χ − χ∂ µ χ * and m * = g * f with g * a coupling. A sum over the fermion generations is understood. The operators in the first line preserve the DM shift symmetries, 5 whereas those in the second line explicitly break them. We assume CP and custodial invariance: the former implies that all the coefficients in Eq. (5) are real, whereas the latter forbids the operator (χ * This Lagrangian is equivalent to the one presented in Refs. [42,43], with a different choice for the operator basis. Our normalization is such that the c and b coefficients have maximal size of O(1), according to the SILH power counting [44]. 6 They can, however, be smaller, for example due to additional symmetries, or because they are suppressed by the degree of compositeness of the SM fermions. Examples of both situations will be discussed momentarily. In Eq. (5) we have neglected additional operators that are expected on general grounds, but do not involve the DM. Among these the operator c H (∂ µ |H| 2 ) 2 /(2f 2 ) is especially important, because it causes a rescaling of all the Higgs couplings by For m χ 10 GeV, the null results of direct searches for DM scattering on nuclei set constraints on the EFT parameters. Only a few directions are strongly constrained. Using the cross section reported in Eq. (21) from Appendix B and taking as example m χ = 100 GeV, the latest XENON1T bound σ χN 10 −46 cm 2 (at 90% CL) [1] translates to where in the first case we set b 1 for the first generation, whereas in the last two cases we took one operator at a time. Naively, the first constraint in Eq. (6) appears to be the strongest. However, we should recall that the expected size of the b coefficients is b Ψ ∼ 2 Ψ , where Ψ is the compositeness fraction of the corresponding fermion field, which is typically very small except for Figure 3: Left: contours obtained requiring that the thermal relic density of χ matches the observed one, assuming the given values of the EFT coefficients and setting the others to zero. Right: the present-day annihilation cross section to SM particles, calculated along the relic density contours shown in the left panel. The black line (yellow band) is the 95% CL observed upper limit (95% CL uncertainty band on the expected limit) from the dSphs analysis in Ref. [46]. The dashed black line corresponds to the observed limit from the analysis of a smaller dSphs sample [47]. The quoted experimental bounds assume annihilation to bb.
the quarks of the third generation: assuming comparable compositeness for the left and right fermion chiralities, , which plugging in the numerical values of m u,d gives f O(100) GeV for any interesting g * . Hence, this bound is easily irrelevant in realistic models. The constraint on m * /c 1/2 B in Eq. (6) (which is due to photon exchange between the DM and the quarks) is more interesting, as it pushes g * toward the fully strongly coupled regime for f ∼ TeV. However, we note that all operators of Eq. (5) that contain the DM current are absent if hidden-charge conjugation, which transforms χ → −χ * but leaves the SM invariant, is preserved. This is the case for all the models presented in Ref. [15]. Furthermore, these operators vanish trivially for real scalar DM. More important is the constraint on λ in Eq. (6), which has crucial implications for model building, even though this coupling only arises at loop level for pNGB DM: models where λ is generated by top loops [2,4,12] have either been excluded or are presently being tested at XENON1T. However, recently Ref. [15] proposed several scenarios where λ is very suppressed, either due to the smallness of the light SM fermion masses or because it arises at higher loop order. In these models the DM-nucleon scattering cross section is well below the current bound, and may even be under the neutrino floor.
While our original motivation for the Lagrangian in Eqs. (4,5) are theories where both the Higgs and the DM arise as composite pNGBs, this EFT can also apply to models such as that of Ref. [17], where a complex scalar S is added to the SM as The U (1) symmetry acting on S is broken spontaneously by the VEV v s , where S = (v s +σ)e iφ/vs / √ 2 , and explicitly by the term proportional to µ 2 S , yielding φ as a real pNGB DM candidate with mass µ S , stabilized by the remnant S → S * invariance. If the radial mode σ is sufficiently heavy, it can be integrated out, obtaining an effective Lagrangian that after a field redefinition matches Eq. (5) with c d /f 2 = λ HS /m 2 σ . 7 In this setup λ arises at one loop and results in a very suppressed cross section for DM-nucleon scattering [18,19]. Gravitational wave signals from phase transitions in the early Universe have also been explored [22,26]. Extended models of this class were discussed in Refs. [21,23,25].
In summary, there exists a class of WIMP models where the direct detection constraints are structurally satisfied due to the pNGB nature of the DM, and in particular λ has negligible impact on the phenomenology. Depending on the specific realization, the DM-nucleus scattering signal may be within reach of future experiments, mediated for example by the contact interactions between the DM and the down quarks parametrized by c χ d , if the DM mass arises from bottom quark loops [15]. However, the signal may also be hidden below the neutrino floor, making its detection extremely challenging. The key element all these scenarios have in common is the derivative Higgs portal operator c d , which is negligible in DM-nucleus scattering, but unsuppressed in annihilation. The observed relic density is reproduced via the freeze-out mechanism for f ∼ TeV, as shown in the left panel of Fig. 3. For completeness we show results down to f = 100 GeV, although it should be kept in mind that for weak-scale f , finding ultraviolet models that are acceptably described at low energies by Eq. (5) may be challenging.
Since c d mediates s-wave annihilation, indirect searches are crucial probes of pNGB DM. In the right panel of Fig. 3 we compare the present-day annihilation cross section to the bounds obtained by the Fermi-LAT [46,47] from gamma-ray observations of dwarf spheroidal galaxies (dSphs). For DM mass of O(100) GeV, these bounds (which were computed assuming DM annihilation to bb) are very close to the canonical thermal cross section σv rel can ≈ 2 · 10 −26 cm 3 s −1 . We see that the region m χ 70 GeV is already in mild tension with data, and future improvements of the sensitivity will provide a stringent test of pNGB DM. Furthermore, searches for DM annihilation to monochromatic photons can be relevant for m h /2 m χ m W , where the branching ratio of the off-shell Higgs to γγ is 10 −3 . We find that in this region the expected cross section is close to the Fermi-LAT sensitivity [48], if a favorable DM profile (contracted Navarro-Frenk-White) and the accordingly-optimized region of the sky are considered.
The recent measurement of the antiproton spectrum by AMS-02 [49] gives additional constraints, as well as intriguing hints of an excess [50,51]. Systematic uncertainties are larger in the antiproton channel, although they have been argued to be under control [52]. Recently, an explanation of the excess as originating from annihilations of pNGB DM has also been proposed [24]. As the situation has not settled yet, we do not show antiproton results here.
Clearly, should a robust excess emerge in indirect detection, its verification by other probes would be paramount to establish its origin. For pNGB DM, collider searches constitute the most important test, given the structural suppression of the direct detection cross section. The key ingredient of the scenario is the derivative Higgs portal, which can be probed in DM pair production as studied in this paper. Additional signatures may arise, depending on the values taken by the b coefficients in Eq. (5). A theoretically motivated situation is one where the b's are large only for the third generation quarks. 8 Then the process gg → tt + / E T becomes relevant; at one loop, a contribution to monojet is also generated. The interplay of these two processes at the LHC was studied in Ref. [53] for fermionic DM. Here we limit ourselves to note that for scalar DM the tt → χχ * amplitude at high energies is not suppressed by an m t insertion, plausibly leading to enhanced sensitivity.

The off-shell Higgs portal in vector boson fusion
We consider φ pair production via off-shell Higgs in VBF, which proceeds through the diagram in Fig. 1. The cross section is written as (V = W or Z) whereŝ = M 2 φφ . The parton luminosity is The parton distribution functions (PDFs) for longitudinal polarizations, which are the only ones coupled to the Higgs, read in the limitŝ where for the W , C f v = − C f a = g/(2 √ 2), and for the Z, For the derivative and marginal Higgs portals, the partonic cross sections arê where we took the high-energy limitŝ m 2 V . Notice the relative enhancement of the derivative portal at high energy,σ deriv /σ marg ∝ŝ 2 . Henceforth we describe our analysis, considering in turn high-energy lepton colliders and hadron colliders.

High-energy lepton colliders
At lepton machines the collision energy is equal to the collider center of mass energy, apart from the effect of initial state radiation, which we neglect in this paper. Assumingŝ m 2 h , we can perform the integral in Eq. (8) and obtain a simple analytical expression for the cross section. For the derivative coupling we find where R eē w ≈ 0. Note the very different scalings with the collider energy and DM mass: for s m 2 φ the derivative portal gives σ ∝ c 2 d s/f 4 , which grows very quickly with √ s and is almost independent of m φ , whereas the marginal coupling leads to σ ∝ λ 2 log(s/m 2 φ )/m 2 φ , weakly dependent on the collider energy and rapidly decreasing as the DM mass is increased. Since the partonic cross section receives an important contribution from the threshold regionŝ ≈ 4m 2 φ , Eqs. (13) and (14) are accurate only if m φ m h /2. In our analysis we always use the exact cross sections as calculated by MadGraph5 [56] using a FeynRules 2.0 [57] implementation of Eqs. (1) and (2).
At lepton colliders we use MadGraph5 to perform a parton-level analysis, which is a reasonable approximation given the extremely clean final state consisting of two leptons and missing momentum. The dominant background is e − e + → ννe − e + . At generation we require p e T > 10 GeV, |η e | < 5 and ∆R ee > 0.4. The most useful kinematic variables to separate signal and background are: the missing transverse energy and missing invariant mass, defined as where / p = ( √ s, 0 )−p e − −p e + ; as well as the invariant mass M ee and pseudorapidity separation ∆η ee = |η e − − η e + | of the electron-positron pair. Distributions for the last three variables are shown in Fig. 4 for CLIC 3, chosen as representative example, and Table 4 describes the flow for the optimized cuts we adopt. For the derivative portal M ee provides less discriminating power compared to the marginal portal, but this is compensated by a tighter requirement on the MIM. After all cuts the background is dominated by ν eνe e − e + , including a contribution of O(10)% from on-shell W W production.   As the derivative Higgs portal is a non-renormalizable operator, we must pay attention to follow a consistent procedure when setting limits on its coefficient. For this purpose, we adopt the method proposed in Ref. [58]: in addition to applying the selection cuts we discard events with characteristic energy E > g * f /c 1/2 d , where g * ≤ 4π is a coupling, because such events cannot be described within the EFT. We choose E = M φφ and for a given value of g * we derive an exclusion region in the plane (m φ , f /c The analysis at a future muon collider is similar, except that at the very high center of mass energies under discussion for such a machine [59] (we take √ s = 6 and 14 TeV as benchmarks, see Table 1), the muons produced by the signal have extremely large pseudorapidity. The problem is most  Only for this plot, the generation-level cut on the muon pseudorapidity was relaxed to |η µ | < 8. Middle: projected 95% CL constraints on the derivative portal at FCC 100 and at a muon collider, assuming the forward detector coverage extends up to |η| = 6 (solid) or |η| = 4.7 [5] for the FCC [muon collider] (dashed). Right: same as the middle panel, for the marginal portal. severe for √ s = 14 TeV, illustrated in the left panel of Fig. 5. If the detector coverage is limited to |η µ | < 5 (corresponding to ∆η µµ < 10) an O(1) fraction of the signal is lost, causing a significant degradation of the sensitivity. Extending the detector capability in the forward region to |η µ | < 6 removes this problem, capturing almost all the signal with a very mild increase in background. The resulting gain in sensitivity is shown in the middle and right panels of Fig. 5. The muon collider projections in Fig. 2 assume the extended coverage |η µ | < 6.

Hadron colliders
To obtain the cross section at hadron colliders we convolute Eq. (8) with the parton luminosities in the proton, with where √ S is the collider energy, Q the factorization scale, and the sum in Eq. (16) extends over all relevant combinations of quarks and antiquarks. For m φ m V /2 this expression provides a sensible estimate of the cross section: we have checked that already for m φ = 100 GeV the agreement with the exact computation is within 30%, improving to < 10% for m φ > 200 GeV.
Our HL-LHC analysis is an extension of the recent searches for VBF-produced, invisibly-decaying Higgs at CMS [38] and ATLAS [62]. We generate events at parton level in MadGraph5 aMC@NLO, shower them with Pythia8 [60] and pass them through fast detector simulation with the help of Delphes3 [61], using the CMS card. The main backgrounds are Z νν +jets and W ν +jets. Each of them is separated into a QCD part, where the jets arise from strong interactions, and an electroweak (EW) part, where the jets are emitted purely via weak couplings. The interference between the two components is negligible [62]. Despite the smaller total cross section, the EW contributions are important because their kinematics closely resembles that of the signal, so that they constitute an O(1) fraction of the total background after all cuts are applied. For the QCD processes we generate V + 2 jets at NLO, whereas for the EW contribution we use a LO matched sample of V + 2, 3 jets. We fix the overall normalization of each sample by comparing our 13 TeV predictions with the expected event yields of the CMS shape analysis, reported in Table 3 of Ref. [38]. After this rescaling the shape of all samples agrees within 10% with the CMS expectation, as shown by the left panel of Fig. 10 in Appendix A. We then apply the same rescaling factors to our 14 TeV samples. We also include the tt+jets background, which is generated at LO matched with up to 2 additional partons, and normalized to the NNLO+NNLL cross section of 974 pb [63].
The baseline cuts we adopt in our 14 TeV analysis are Normalized distributions for the signal and background after these initial cuts are shown for m φ = 100 GeV in Fig. 6. An important discriminating variable is ∆η jj , which is most effective for the derivative portal (Fig. 6 -left). This is the opposite of what happens at lepton colliders, where the ∆η distribution is harder for the marginal portal, see e.g. the left panel of Fig. 5. We fix ∆η jj > 4.8 [4.2] for the derivative [marginal] portal. Additionally, we consider a more stringent cut on / E T and a cut on m jj . The significance obtained by varying these two requirements is shown in Fig. 7. It is apparent that a relatively mild cut on / E T would be preferred for both portals. However, as events must be selected using a missing energy trigger, a rather stringent requirement on / E T needs to be imposed. Following the ATLAS analysis [62] we set / E T > 180 GeV. Note that the shape of the missing energy distribution is very similar for the signals and the EW backgrounds (middle panel of Fig. 6). Finally, we optimize the m jj cut, separately for the derivative and marginal signals. The complete cut flow for the benchmark m φ = 100 GeV is reported in Table 5. The tt+jets background is small, and we   neglect it in our HL-LHC projections. Note that the EW V +jets are very important, making up 40% of the total background to the derivative signal. For different m φ hypotheses all cuts are kept identical, except for m jj , which is optimized in each case. It is interesting to check the effect on the significance of varying the missing energy cut. We take as examples / E T > 150 GeV, which at ATLAS corresponds to a trigger efficiency of ∼ 90% on signal events [64] (to be compared with 98% for the reference choice / E T > 180 GeV [62]) and / E T > 250 GeV, as required in the CMS analysis [38]. For the derivative portal with m φ = 100 GeV the bound on f /c 1/2 d varies from 290 to 260 GeV, whereas for the marginal portal with naturalness-inspired coupling strength the reach on m φ varies from 135 to 117 GeV.
Our FCC analysis proceeds along similar lines. We apply to the V +jets backgrounds the same rescaling factors derived from the comparison with CMS at 13 TeV, whereas tt+jets is normalized to the theoretical prediction of 34.7 nb [65]. For Delphes3 we use the FCC card. The baseline cuts are the same as in Eq. (18), except for the requirement on the jet pseudorapidity, which is increased to |η j 1 ,j 2 | < 6. As emphasized in the FCC-hh Conceptual Design Report [66], extending detector coverage up to |η j | = 6 would be important for the measurement of VBF processes. This is quantified for our signal in the middle and right panels of Fig. 5, where we show the gain in sensitivity obtained extending the forward jet measurement from the LHC range |η j | < 4.7 to the expected FCC design |η j | < 6 (note that the FCC 100 bounds shown in Fig. 2 assume coverage up to |η j | = 6). Additionally, we set ∆η jj > 5.5 [4.2] for the derivative [marginal] portal, whereas the missing energy cut is fixed to / E T > 200 GeV [41]. The complete cut flow for m φ = 100 GeV is given in Table 6. For the marginal portal the tt+jets background is not negligible, hence we consistently include it in our FCC projections.  For the HE-LHC we again apply the 13 TeV rescaling factors to the V +jets backgrounds, whereas tt+jets is normalized to 3.73 nb [67]. We use the HL-LHC Delphes card. The baseline cuts are as in Eq. (18), except that we impose |η j 1 ,j 2 | < 4.9. As additional requirements we apply ∆η jj > 5.2 [4.2] for the derivative [marginal] signal, fix / E T > 200 GeV, and as before we optimize the m jj cut for each m φ and portal hypothesis.
A comment is in order about the role of systematic uncertainties, which are consistently neglected in our analysis, since our aim is to show the best sensitivity that can be achieved in principle. At lepton colliders the S/B ratio is at least several per-cent (see e.g. Table 4), therefore the effect of systematics is expected to be minor. By contrast, at hadron colliders the S/B is below the per-mille, as can be read in Tables 5 and 6, hence systematics will have an important impact on our hadron collider limits. Its assessment requires a careful analysis that should be done at the experimental level, taking advantage of the theoretical advancement in the precise predictions of the dominant V +jets backgrounds [68].

Discussion and outlook
We begin our discussion with a few comments about Fig. 2. An important result is that hadron colliders have stronger sensitivity to the derivative Higgs portal than to the marginal Higgs portal, because the former produces harder kinematic distributions that allow a more effective background suppression (see also Tables 5 and 6). To make this statement more quantitative, we show in Fig. 8 the bounds on the (absolute value of the) effective coupling strength at the kinematic threshold M φφ = 2m φ , which is c d 4m 2 φ /f 2 and λ, respectively. For hadron colliders the threshold is a reasonable point of comparison, due to the PDF suppression at higher M φφ . We find a better reach on the derivative portal, thanks to the tail at higher invariant masses. For m φ = 100 GeV, the ratio of the effective coupling bounds is ∼ 1/4.0 at the HL-LHC, ∼ 1/5.5 at the HE-LHC and ∼ 1/11 at the FCC. A consequence of this increased sensitivity is that, for example, FCC 100 gives projected constraints comparable to a 6 TeV muon collider for the derivative coupling, but very similar to CLIC 3 for the marginal coupling.
Our key motivation to consider the derivative portal is pNGB DM. To assess the future prospects we follow the relic density contours in the left panel of Fig. 3, where we read that m χ > m h /2 corresponds to f 500 GeV (for complex DM). This region is out of the reach of HL-LHC and CLIC 1.5, can be probed to a very limited extent at the HE-LHC and CLIC 3, and would be truly explored only at FCC 100 or a muon collider. In particular, a µC 14 would test DM masses up to about 600 GeV.
As a representative example of the marginal portal, we consider scalar top partners with coupling fixed by naturalness. In this case even FCC 100 will have a limited sensitivity, probing top partner masses up to approximately 300 GeV. A muon collider would have a clearly superior reach, extending all the way up to mũc ≈ 1 TeV for center of mass energy of 14 TeV. This would constitute an impressive, model-independent test of the possibility that the little hierarchy is stabilized by scalar top partners, regardless of their decay pattern.
We conclude with some comments about other potential probes of the pNGB DM scenario. Firstly, in this setup the couplings of the Higgs to visible particles are generically expected to show deviations from the SM predictions. The size and pattern of the corrections is, however, dependent on the model (and in particular, on whether the Higgs also arises as a pNGB or not), and their discussion is outside the scope of this paper. As far as probes of the pNGB DM itself are concerned, we emphasize that the strategy followed in this first study is simply the most direct one, exploiting the tree-level production of a DM pair through an off-shell Higgs. Probes of virtual effects may provide competitive sensitivity, in analogy to the one-loop tests of the marginal Higgs portals proposed in Refs. [69][70][71]. We believe the study of such indirect sensitivity to pNGB DM deserves further attention [72].
Acknowledgments We have benefited from conversations with R. Frederix, J. Herrero-García and P. Panci. We also thank S. Hong for correspondence about Ref. [32], and A. Tesi for clarifications about Ref. [55]. This research has been partially supported by the DFG Cluster of Excellence 153 "Origin and Structure of the Universe" and by the Collaborative Research Center SFB1258. MR is supported by the Studienstiftung des deutschen Volkes, and acknowledges the hospitality of the Cornell theory group in the final stages of the project. All authors warmly thank the MIAPP, where part of this work was done.

A Details of the analysis and comparison with previous works
At lepton colliders we cut on four kinematic variables, namely MIM, ∆ , / E T and M : • The cuts on the MIM are summarized in the left-most panel of Fig. 9 for the derivative portal.
For the marginal portal we always take MIM > 2m φ , and in a few cases we also impose an upper   • The cuts on ∆η are shown in the middle-left panel of Fig. 9. For the marginal portal at ILC 1 this cut is replaced by H T (ee) < 260 GeV, as in Ref. [32].
• The cuts on / E T = MET are shown in the middle-right panel of Fig. 9 (the points that are not immediately visible in this plot all have / E T > 80 GeV, e.g. the marginal portal at CLIC 3).
• The cuts on M are shown in the right-most panel of Fig. 9.
For hadron colliders, we show in the middle panel of Fig. 10 the cuts on m jj .

A.1 Comparison to Ref. [36]
We deem it useful to provide a thorough comparison of our results for the marginal portal at hadron colliders to those of Ref. [36]. For this purpose we have tried to reproduce their VBF analysis as closely as possible, given the available information. The results are shown in the right panel of Fig. 10, together with the bounds quoted in Ref. [36] and those obtained in our main analysis. 10 We are able to reproduce well the results of Ref. [36], agreement being excellent for the FCC and good for the HL-LHC. This exercise helps us to pinpoint the differences between our analysis and that of Ref. [36]: 1) In Ref. [36] the EW components of V +jet were neglected, while the QCD components were generated at LO and normalized to the MadGraph5 cross sections. In this work we have included EW production and generated the QCD part at NLO; normalizations were fixed by comparison with the 13 TeV CMS results [38]. 10 Only for this plot, following Ref. [36] we adopt the coupling normalization c φ = λ/2 and use a two-sided test to determine the exclusion bounds, i.e. we require S = √ 2 erf −1 (1 − p) ≈ 1.96 for p = 0.05 (95% CL).
2) In Ref. [36] both the / E T and m jj cuts were optimized, without imposing a trigger-motivated lower limit on the former. As a result the optimal values were / E T 120 GeV at the HL-LHC and / E T 100 -160 GeV at the FCC. In this work we have fixed / E T > 180 and 200 GeV, respectively, and optimized only the m jj cut.
3) In our FCC analysis we have assumed forward jet tagging up to |η| = 6, to be compared with 4.7 used in Ref. [36].
Owing to the points 1) and 2) our HL-LHC sensitivity is weaker, while at the FCC the point 3) more than compensates for the first two, allowing us to derive slightly more stringent limits.

B Cross sections for pNGB Dark Matter
To calculate the cross section for DM scattering on nuclei, we first match Eq. (5) to the following effective DM-SM interactions, ψ=u,d,c,s,t,b m ψ a ψ |χ| 2 ψψ + iλ V ψ (χ * ↔ ∂ µ χ) ψγ µ ψ + iλ A ψ (χ * ↔ ∂ µ χ) ψγ µ γ 5 ψ + ie m 2 * c B (χ * ↔ ∂ µ χ)∂ ν F µν (19) with coefficients From these we derive the spin-independent DM-nucleon cross section, where δλ V ψ = c B e 2 Q ψ /m 2 * and m N = (m p + m n )/2. For the scalar operators, we have separated the contribution of the light quarks from that of the heavy quarks. The associated form factors take the values (see App. D of Ref. [15] for the relevant references) f p Tu = 0.021, f p T d = 0.041, f n Tu = 0.019, f n T d = 0.045, and f p,n Ts = 0.043, leading to f p,n Tg = 1 − ψ=u,d,s f p,n T ψ 0.89. It is important to note that for the vector current operators proportional to λ V ψ , the form factors simply count the valence content of the nucleons, namely f p Vu = f n V d = 2 and f p V d = f n Vu = 1. The additional term proportional to δλ V ψ originates from t-channel photon exchange. As for DM annihilation, the last two operators in the first line of Eq. (5) contribute to the p-wave, while all other operators contribute to the s-wave. For m χ m W the χχ * → W W, ZZ, hh channels dominate and the b Ψ and c B operators are subleading to the unsuppressed Higgs portal. Conversely, for m χ m W the χχ * → ψψ channel (where ψ is a SM fermion) dominates and the b Ψ and c B operators are important at freeze-out, when their velocity suppression can be less severe than the m ψ -suppression that characterizes the other operators. As an example we consider χχ * → bb, for which we find a cross section where c hbb = g hbb /g SM hbb and the coupling of the SM fermion ψ to the Z is defined as ig Z γ µ (v ψ − a ψ γ 5 ). Taking the thermal average and considering for simplicity the limit m χ m b we arrive at σv rel (T ) N c 4π m 2 b A 2 (4m 2 χ ) + 4m χ T B 2 (4m 2 χ ) + C 2 (4m 2 χ ) .