Light Dark Sectors through the Fermion Portal

Pairs of Standard Model fermions form dimension-3 singlet operators that can couple to new dark sector states. This"fermion portal"is to be contrasted with the lower-dimensional Higgs, vector and neutrino singlet portals. We characterise its distinct phenomenology and place effective field theory bounds on this framework, focusing on the case of fermion portals to a pair of light dark sector fermions. We obtain current and projected limits on the dimension-6 effective operator scale from a variety of meson decay experiments, missing energy and long-lived particle searches at colliders, as well as astrophysical and cosmological bounds. The DarkEFT public code is made available for recasting these limits, which we illustrate with various examples including an integrated-out heavy dark photon.


Introduction
In recent years, the search for dark matter has broadened beyond weakly-interacting massive particles in the GeV-TeV range. There has been a resurgence of interest in improving experimental sensitivity to sub-GeV dark matter, while much progress has also been made in neutrino experiments (see e.g. Refs. [1][2][3] for recent overviews). In addition to uncovering the nature of dark matter, these searches could open a window onto a rich dark sector that must often accompany it. The dark sector's experimental signatures often share many similarities with those of dark matter and neutrinos. Moreover, dedicated experiments have been proposed to look specifically for the spectacular signal of longlived particles decaying back to the visible sector. The generality of such dark sectors requires an appropriate framework in which to characterise the sensitivity of these various searches. Many models, either simplified or top-down motivated, have been constructed 1 and used as benchmarks. Here we propose a bottom-up effective field theory approach to characterise the phenomenology of light dark sector searches more generally. Dark sector fields are singlets under the Standard Model SU (3) c × SU (2) L × U (1) Y gauge groups. They are motivated by observations of neutrino masses and dark matter, which require additional particle content that may well involve an extended dark sector, and arise generically in many models of new physics addressing a variety of problems. Additionally, dark sectors with potentially rich phenomenology are a generic prediction of compactified string theory (see e.g. [4][5][6][7] for recent studies and reviews). Indeed, there is nothing exotic about singlet charge assignments under the Standard Model gauge groups; we already know of particles that are SU (3) c × SU (2) L singlets, so it is reasonable to suppose others may exist that go a step further in being uncharged under U (1) Y as well.
In full generality, the visible sector of the Standard Model can be described as communicating with dark sectors through so-called "portal" operators, O (1.1) The quantity Λ ij is a dimensionful scale and c ij is a dimensionless coefficient; if d + d > 4 then the interaction is associated with an effective non-renormalisable Lagrangian where Λ ij is related to the heavy mediator mass and c ij is a Wilson coefficient. If the light dark sector is only connected to the Standard Model through heavy mediators, as for example in hidden valley models [8], the resulting scale suppression would provide a natural explanation for the weakness of the dark sector's interactions with the visible sector.
The first few portals ordered by dimensionality are listed in Table 1. The lowestdimensional portal operators are the well-known Higgs, vector, and neutrino portals (see e.g. Refs. [9][10][11][12][13] and references therein): where d is the mass dimension, F µν is the electromagnetic (or hypercharge) gauge field strength and the Higgs and lepton doublet fields are denoted by H and L. The low operator dimensionality of these portals allows them to form renormalisable interactions (d + d ≤ 4) with the hidden sector. For example, the Higgs portal may couple to another scalar; the vector portal can kinetically mix with a hidden sector gauge field strength; and the neutrino portal can be responsible for neutrino masses through a right-handed neutrino coupling. The next lowest-dimensional singlet operators involving Standard Model fermions are those that we dub "fermion portals". These are singlet combinations of Standard Model fermion fields, of which the lowest-dimensional takes the form where ψ represent the Standard Model fermions and the different types of contractions are left implicit. At the level of unbroken electroweak symmetry, these portals can only form higher-dimensional operators of dimensions 5, 6 or greater with the dark sector, 1 suppressed by some effective field theory cut-off scale. The heavier the cut-off, the weaker the interaction between the visible and dark sectors. Searches for light, weakly-coupled dark sectors in fermion portals could then also yield complementary information about the scale of heavier new physics mediating the interaction. For example, the neutral current dimension-3 fermion portal of Eq. 1.5, O SM , can form a dimension-5 operator with a derivatively coupled scalar, SM . (1.6) The prototypical example of this is the familiar axion, where the scale suppression is related to the axion's symmetry-breaking scale. For some recent phenomenological studies of constraints on the scale of the axion decay constant, see for example Refs. [14][15][16][17][18].
In this paper we characterise the sensitivity of dark sector searches focusing on the case of neutral current fermion portals to a pair of light dark sector fermions. 2 In this case the dark sector fermion χ forms a dimension-6 four-fermion operator with the Standard Model fermion pair of Eq. 1.5, where we focus on the tensor structures Γ ⊂ {γ µ , γ µ γ 5 }. We consider in this work fourfermion operators coupling the dark sector fermions to quarks (with effective couplings g ij u and g ij d where the indices i, j refer to the generation of SM fermions) and charged leptons (with effective couplings g ij l ). We will typically choose particular ratios g u : g d : g l and present the limits in terms of Λ/ √ g, where Λ is a mass scale and g represents the overall coefficient (that we also refer to as the effective coupling). 1 The portal of Eq. 1.5 can make a dimension 4 interaction with a hidden sector vector boson, but a light gauge boson is conventionally included as part of the vector portal phenomenology. Here we shall be concerned with heavy mediators. 2 We note that there are also 3-fermion singlet combinations, ψ i ψ j ψ k (d = 9/2), such as d i u j d j . Since they carry baryon number B = 1 they must couple to a hidden sector fermion that carries the opposite baryon number. Such a singlet fermion with baryon number and no lepton number has been discussed e.g. in Refs. [19,20]. Ref. [19] also categorises the various singlet fermion lepton number assignments that would forbid the renormalisable neutrino portal operator at tree level. Indeed, there exists a rich set of possibilities for dark sector fermions beyond the familiar right-handed neutrino with lepton number L = 1, though we will not consider them any further here.
Light dark sectors are typically probed by an extremely wide range of experiments usually referred to as the intensity frontier of particle physics. They share a relatively low center-of-mass energy, high available statistics and very good background rejection. We will use three general type of accelerator-based experiments (a complete list of the searches we have implemented can be found in Table 7 in Appendix B): • First, there are the dedicated flavour experiments, which can be either based at e + e − colliders such as the B-factory experiments BaBar [21], Belle-II [22], or at beam dump facilities such as the Kaon factories NA62 [23] or E949 [24]. These experiments are typically used for missing energy (mono-photon) dark sector searches or for indirect limits using invisible B/K meson decays.
• Second, we have neutrino-focused experiments that are typically based on proton beam dumps. These include the past experiments LSND [25], CHARM [26] or the current near-detector experiments MiniBooNE [27], and NOνA [28]. Here dark sector particles, abundantly produced at the beam dump, can travel alongside neutrinos and either scatter or decay in the detector if they are sufficiently long-lived.
• The third class of experiments are dark sector-oriented ones. Past experiments were typically searching for axion-like particles or dark photons and were based on electron beam dump experiments. Their sensitivity relies on the electron's bremsstrahlung into the dark sector, whose large cross-section typically compensated for their somewhat lower statistics. Additionally, a significant fraction of proposed new experiments will be either LHC-based (such as FASER [29] and MATHUSLA [30]) or based on a proton beam dump (such as SHiP [31] or a possible extension of SeaQuest [32]). While we aim for a comprehensive coverage of existing experimental limits, we do not attempt to provide projections for all upcoming intensity frontier experiments and will instead focus on a few representative examples. A more complete list can be found in, e.g. Refs. [1,2].
We also discuss astrophysical and cosmological constraints on dark sector fields linked to the SM by a fermion portal. While such particles can constitute all or a fraction of the observed dark matter relic density and must not overclose the universe, 3 we will not derive limits based on this criteria here, as they depend strongly on the inner dynamics of the dark sector states. Indeed, we do not require that the dark sector states make up any of the dark matter, or be cosmologically stable. On the other hand limits from the observed cooling rates of supernovas and the cosmic microwave background can provide relatively model-independent constraints on fermion portal dark sectors with masses below the tens of MeV range. These indirect observational constraints can be complementary dark sector probes to the direct experimental searches listed above.
The paper is structured as follows. In Section 2, we summarise the various production mechanisms relevant for the fermion portal operators, from light to heavy meson decay processes and direct parton-level production. Section 3 focuses then on the various relevant search channels at accelerator and beam-dump based experiments, including decay and scattering signatures of dark sector particles, invisible meson decay limits and mono-X searches. Section 4 presents a selection of the relevant astrophysical limits, from SN1987A cooling to an estimate of early universe cosmological limits. Finally, Section 5 is dedicated to a presentation of the results obtained using our public code DarkEFT. We show different representative choices of models and in particular illustrate the effectiveness of the formalism through the example of a relatively heavy dark photon (in the tens of GeV range). The appendices contain more detail about the meson decay production mechanisms as well as a brief presentation of the DarkEFT companion code, released alongside this paper and dedicated to the recasting of existing dark sector limits into constraints on the fermion portal operators.
2 Dark sector production through the fermion portal Dark sector fermions can be produced via the fermion portal through various processes involving quark and lepton bilinears. In this work we will be agnostic as to the flavour pattern of the Wilson coefficients of the effective operators. These may in principle involve both flavour-diagonal and off-diagonal couplings. Our results will be as general as possible, so that they may be applied to any model of dark sector fermions coupling to the SM via a mediator far off-shell.
The production channels of dark fermions we consider here are: • Light meson decays: this channel depends on the nature of the meson and the operator -axial vector or vector. It typically proceeds through an associated decay, e.g. π 0 → γχχ, or via a fully dark decay, e.g. π 0 → χχ. This channel requires non-zero couplings g u or g d to up or down quarks respectively.
• Heavy meson decays: this channel will be important for vector charmed quarks, e.g. J/Ψ, or for flavour-violating rare meson decays. This channel depends on couplings to heavy quarks such as g c to charm or g b to bottom. The latter in particular will require non-flavour-diagonal couplings.
• Direct production either at the parton level, via pp → χχ or pp → jet + χχ, or in electron colliders, via ee → χχ or mono-photon ee → γχχ.
The production of light dark sector fermions through bremsstrahlung, p(e)N → p(e)N χχ, is also possible, although often with smaller cross-sections than some of the processes listed above. This production mechanism, when paired with a search for missing momentum/energy as an experimental signature, can be quite powerful to search for light dark sector particles [37][38][39], as it does not require a visible decay, and therefore a further g/Λ 2 suppression. We defer study of this production mechanism and detection approach to future work.
Interestingly, parton-level production is usually irrelevant in standard portal searches for low mass dark sectors, since QCD perturbativity breaks down before the mediator scale. In our fermion portal model, this is no longer the case as long as the mediator mass is larger than the QCD scale, so that direct production becomes essentially constant for low dark sector masses.
In the following, we will consider the most generic case in which there are two dark sector states, χ 1 and χ 2 , with masses M 1 and M 2 respectively (we will generally take M 2 > M 1 ). All our results can be trivially applied to a one state scenario by setting M 1 = M 2 . We now consider in turn each of these production mechanisms.

Light meson decays
Light mesons are abundantly produced in proton-based colliders and beam dump experiments. Furthermore, their decays typically proceed with a relatively long lifetime which tempers the strong suppression from the fermion portal's high scale. For a meson M , the final production number N prod of the dark sector state χ 2 (in association with χ 1 and possibly another SM particle X) is given by where M m is the meson mass and Λ is the scale of the fermion portal operator. Heavy mesons are then expected to interact more with the dark sector than lighter fermions and can dominate the production rate. This is somewhat reminiscent of Higgs portal phenomenology, where the dark sector also couples more strongly to heavier fermions, though the effect is even more pronounced in the fermion portal case as it has a quartic dependence on the meson mass compared to a quadratic dependence for the Higgs portal. Note that this effect is balanced in part by the fact that lighter mesons have a smaller decay width, as summarised in Table 2, so that their branching ratio into dark sector fields is enhanced. Depending on the pseudo-scalar or vector nature of the light unflavoured mesons, the type of operator (vector or axial vector) will be critical in determining the possible channels for dark sector production in meson decay. We shall consider effective vectorvector (V-V) operators of the form and axial-vector couplings (AV-AV) corresponding to operators of the form where we have included the possibility of having two states in the dark sector, with a mass splitting ∆ χ ≡ |M 2 | − |M 1 |. Depending on which one is the most relevant, we will also use the normalised splitting The splitting can be taken to zero to recover the single state case. We will consider Dirac fermion dark sector states χ 1 and χ 2 . In the following we shall summarise the results for the various amplitudes; details of the calculation can be found in Appendix A. For most decays, we present both the small splitting limit, δ χ 1, and a saturated splitting where M 1 M 2 (δ χ 1). For all numerical results we use the full amplitude which are straightforwardly derived once the effective couplings relevant to the corresponding decaying meson are known. In all the rest of this section, and unless explicitly specified otherwise, one can obtain the decay rates for mixed operators AV-V or V-AV by replacing M 1 by −M 1 in the V-V or AV-AV expressions respectively. The result from the V-V operator can also be used for the so-called "pseudo-Dirac case" [40] when χ 1 and χ 2 are Majorana states originating from a single Dirac field. 4 Interestingly, most of the phenomenology (for instance the type of meson decays, the flavour-violation effects, etc.) depends only on the Standard Model part of the operators so that the following results can also give a rough estimate of the constraints one could expect for other dark sectors not considered here, for instance with scalar fields instead of fermions.
Vector coupling Due to the vector nature of the effective operator, the decay of light pseudo-scalar mesons have to proceed through the axial anomaly with an associated photon production. The dominant production mechanism is then π 0 , η, η → γχχ with a decay amplitude of the form and we have used the effective couplings g 2 P defined in Table 2, and taken P ≡ π 0 , η, η and the pion decay constant f π = 130.7 MeV. The strong numerical suppression factor arises in part from the loop-induced axial anomaly and in part from the phase space suppression for this 3-body decay. Furthermore, we recover as expected the scale suppression by M 4 P /Λ 4 . Both effects combined imply that the dominant production mechanism in most cases will in fact be vector meson decays. Indeed, for dark sectors coupling to the SM through a vector current, vector mesons can decay directly into dark sector particles, e.g. ρ, ω → χχ, 4 The inclusion of a small Majorana mass term triggers the splitting of the Dirac fermion into two Majorana fermions. An important subtlety is the fact that if one insists on keeping positive masses from both Majorana fermions, the mixing matrices become complex. In the case of an initial vector coupling for the Dirac fermion coupling, the interactions term then contains both a leading vectorχ 1 γ µ χ 2 interaction and a sub-leading axial-vector oneχ 1 γ µ γ 5 χ 2 . For larger splitting, both interactions become relevant and it is preferable to use instead a negative mass M 1 with a purely axial-vector interactionχ 1 γ µ γ 5 χ 2 . Note that this limit also extends naturally to the case where the Majorana component dominates, in which case both the masses are positive.

Vector current
Axial-vector current Dark photon Γ M (GeV) with a decay width given by where U ≡ ρ, ω, the effective couplings g 2 U are defined in Table 2. While the Λ suppression remains, the two-body nature of the decay and the absence of α em suppression strongly enhance this decay compared to the previous one. Altogether, as can be seen in Fig. 1a and Fig. 1c, the production for the vector coupling case is strongly dominated by the decay of vector mesons.
In Fig. 1a and Fig. 1c we summarise the corresponding branching ratios for the vector case as a function of the dark sector fermion mass. The couplings are chosen to be either aligned to the electromagnetic ones, g u = 2/3, g d = g s = −1/3, for Fig. 1a, or following a "baryonic" coupling g u = g d = g s = 1/3 for Fig. 1c. We have set the splitting at δ χ = 0 and 20 for the solid and dashed lines respectively. In particular, following the scaling presented in Table 2, we see that the production from ρ meson is strongly suppressed in the baryonic regime.
Axial Vector coupling Here the dominant contributions from meson decay arises from the direct decay π 0 , η, η → χχ, with a decay width given by where we have used the effective couplingsg 2 P defined in Table 2 and taken P ≡ π 0 , η, η . Notice that, similarly to the standard calculation of the decay π 0 → νν, the decay ampli-    tude depends quadratically on the dark sector mass M 1 due to the helicity suppression of the decay amplitudes [41].
In Fig. 1b and Fig. 1d we summarise the corresponding branching ratios for the axialvector case as a function of the dark sector fermion mass with a splitting δ χ = 0 (20) in solid (dashed) lines. The couplings are chosen to be either Z-aligned, with g u = 1/2, g d = g s = −1/2, or uniform across the quarks as in the baryonic case. Notice that this latter case corresponds to a pion-phobic regime and strongly suppresses the production rate at small masses.
Finally, one can obtain the overall number of produced dark sector particles at a given beam dump point by factoring in the typical ratio of mesons per Proton-on-Target (PoT) for the various beam dumps and beam energies. These ratios are listed in Table 3. We have checked using the code EPOS-LHC [42] as distributed whithin the package CRMC [43] 9 Experiment HL-LHC (barn) [29] 14 TeV pp L = 3 ab −1 4.3 0.5 0.05 0.5 0.5 Table 3: Beam and target characteristics for various experiments, along with the total number of protons on target (PoT) and the average number of a given meson per proton. E beam is the beam kinetic energy. The SHiP design is not final. References point either to an analysis paper in the case of existing constraints, or to prospective bounds in the case of future experiments. The ratios quoted for NOνA account only for primary mesons (see [44]). For MiniBooNE, we use the ratios and normalisation from the recent off-target analysis [27].
that the ratios between the different types of mesons were accurate, and further completed the table when only partial information was available (e.g. only the number of π 0 ). The overall normalisation (typically given by the number of π 0 per PoT) tends however to vary strongly from analysis to analysis. In the best cases, GEANT4 simulations of the full hadronic and electromagnetic cascade, supplemented by experimental data are used (as for MiniBooNE in [27]). Several other studies use either PYTHIA8 simulations of a pp process with similar center-of-mass energy, or experimental data from pp collision to extract the meson multiplicity, which tend to underestimate the production (hence leading to conservative limits). Intermediate approaches, such as those presented above using EPOS-LHC, simulate the pN process, but do not include the subsequent showers. In our case, we typically do not choose between the various methods since the overall normalisation depends on the analysis that we will use later for recasting limits; certain studies choose only to keep primary mesons while others use the full hadronic shower components. In particular, this is the case for the projection for NOνA from Ref. [44] which chooses to keep only one π 0 meson per proton on target. For the case of SHiP, where we provide estimates for a 10 events reach, we use the overall normalisation to be around ∼ 10π 0 /PoT, based on the EPOS-LHC results (though we note that a PYTHIA8 approach in Ref. [45] found around ∼ 6π 0 /PoT). from the ones of the lighter mesons, the final experimental efficiencies for dark sector particles originating from their decay cannot be inferred from existing decay and scattering searches which do not include them. Consequently, we will follow two complementary directions: We calculate directly the limits from their invisible decays, for which one need not assume any particular detection efficiencies, and we present naive 10-event projections at SHiP as an order-of-magnitude estimate for the reach of searches based on heavy meson decays.
We shall focus on the decays mediated by the vector effective operator, which can be separated into flavour diagonal ones, φ, J/Ψ, Υ → χχ, and flavour-violating ones such as B → Kχχ, πχχ, etc. The expression for the decay width of the former has been already estimated in Eq. (2.6), with the corresponding decay constants given in Table 6 of Appendix A.
For the three-body decay of heavy pseudo-scalar mesons, we extend the study of Ref. [17] for the case of axions to our four-fermion portal operator generated by an offshell mediator (see also Ref. [47]). We consider the B and K mesons, whose sensitivities are enhanced by their heavier mass and flavour-violating decays through off-diagonal vector couplings: The three-body decay of a heavy meson P (mass M ) into a lighter meson P (mass M ) and two dark sector fermions χ through a contact interaction is derived in Appendix A.
In the massless χ limit, denoting g P P as the relevant effective coupling, the decay width We leave for future work the recasting of existing studies such as [46].

11
reduces to The hadronic form factor f + (q 2 ) is obtained from lattice results in Ref. [48]; for a momentum transfer q 2 ≡ (p P − p P ) 2 → 0 its value is between 0.23 and 1 depending on the decay process (see Table 6 in Appendix A). Since the q 2 dependence of the form factor does not vary significantly for the purpose of setting limits, we take these constant values as a good approximation. The branching ratios are plotted for various B and K meson decays in Fig. 2 as a function of the light dark sector mass with no splitting between the two dark sector masses.
In order to get an estimate for the number of B mesons produced at SHiP, we multiply the number of protons on target, N prot = 2 × 10 20 , by the ratio of the B meson production cross-section per nucleon, σ B 3.6 nb, to the total proton-nucleon cross-section, σ pN 40 mb. For K 0 S , K − and K + we find that the multiplicities of 0.232, 0.224 and 0.331, respectively, from Ref. [49] agree well with our EPOS simulations for proton-proton collisions at √ s = 27.4 GeV. However, the target for SHiP is tungsten where we find instead a factor of ∼ 2.5 enhancement in multiplicities. These numbers are used to generate a naive projection of limits assuming 10 events, based on the routine presented in the next sections. 6

Parton-level dark sector production
When the centre-of-mass energy E cm of the relevant beam-dump or collider experiment is larger than O(1) GeV, direct parton-level production can become relevant. In the case of dark sector production at beam dump experiments for a target of atomic number Z and mass number A, the final direct cross-section production is given by The total number of produced pairs of dark sector particles can then be deduced as 11) where N PoT is the number of protons on target and σ p+A is the total scattering crosssection on the material. The second equality makes the (strong) assumption that "screening" effects, which make the typical proton-nuclei cross-sections scale proportionally to A m rather than A (typically with m ∼ 0.7 for the energy of intensity frontier experiments), apply similarly to new physics processes as to the Standard Model 7 (see e.g. [48,51,52]).  Figure 3: (a) Cross-sections for the proton partonic production processes at Λ = 1 TeV for uū, dd → χχ for the SPS (400 GeV) beam and for the NuMI (120 GeV) beam (b) Cross-section times Λ 4 at LHC for pp → χχ + jets as a function of the new physics scale Λ.
Assuming for now that both final dark sector states have similar mass, the crosssection can be written for a process with a center-of-mass energy √ s as where we have introduced the parton distribution function (PDF) f N q for the parton N (proton p or neutron n), the sum runs over all quarks and antiquarks and the initial momentum of the quarks in the cross-section σ(qq → χχ) depends on the momentum fractions x 1 and x 2 .
In practice the neutron PDFs can be determined from those for protons using isospin symmetry. Factorising the effective couplings between the dark sector fermions and the quarks, g u and g s , we obtain where we neglected the quark masses and factored out the couplings in the last crosssection to obtain 8  Table 4: Cross-section for direct production in various high-energy proton beam experiments for M 1 +M 2 E cm . The effective coupling were chosen electromagnetically-aligned and Λ = 1 TeV.
It turns out that in the regime of interest, the PDF ratios follow the scaling where σ 0 p+p→χχ is estimated in the electromagnetic alignment g u = 2/3, g d = −1/3 and shown in Table 4 in the limit M 1 , M 2 E cm . This scaling is accurate at 20% in the region of low M 1 + M 2 mass compared to the center-of-mass energy E cm . We have estimated the cross-section using CTEQ6.6M pdfs both directly from the above formula and through MadGraph5 aMC@NLO platform [53] with the effective theory model implemented in FeynRules [54] to create the UFO module [55]. The renormalisation scale choice was left dynamical, chosen as the sum of the transverse mass of the outgoing dark sector fields divided by 2. 9 We show in Fig. 3a the cross-sections for the uū, dd → χχ (for a proton-proton process), which corresponds to choosing g u = 1, g d = 0 and g d = 1, g u = 0 respectively in the equations above. Note that in the numerical process, we combine both curves using the relations from Eq. (2.17) and account for a slight dependence of the coefficient on the energy of the initial beam.
We implemented the dark photon direct production with the same approach in order to recast the limits later on. In this case the renormalisation scale is set to the dark photon mass corresponding to the resonant Drell-Yan production [12]. Since the dark photon is typically quite light, we probe a different renormalisation scale range in this case, but one still has typically As long as the total mass of the dark sector fermion pair M 1 + M 2 E cm Λ, the cross-section is roughly constant and independent of the actual dark sector masses. But an additional difficulty arises when the effective scale becomes of the order of the centreof-mass energy of the process considered; the effective theory starts becoming unreliable 15) and the couplings are given as g 2 q = g D eQ q with Q q the quark electric charge, g D the dark gauge coupling and the kinetic mixing. 9 When estimated directly, we chose Q = 500 GeV for the LHC though we note that the ratio was not significantly modified by varying it from 250 GeV to 1 TeV. since direct mediator production should dominate. This issue is mostly relevant for LHC searches and has been intensely studied in recent years following searches for dark matter at the LHC through effective operators. While various approaches have been suggested (see e.g. [56] for a brief summary), the general strategy is to restrict on an event-by-event basis at the Monte-Carlo generator level the typical energy of the process to be below the effective field theory scale. 10 While we will later consider briefly for completeness some of the mono-X limits from the LHC, our dominant interest and strongest bound will come from direct "dark" production of a dark sector pair that is subsequently detected by a dedicated experiment such as FASER or MATHUSLA. Removing the requirement for an additional high-p T particle significantly increases the potential reach of the EFT. We illustrate this in Fig. 3b for the associated pp → χχX dark sector production (relevant for mono-X searches) where we present the typical cross-section as a function of the effective operator scale, following the procedure described in Ref. [56].
3 Hunting for the fermion portal's dark sectors 3

.1 Long-lived dark sectors
Searches for hidden particles in both beam-dump and accelerator-based experiments can typically be decomposed between a production stage and detection stage. The former was described in the previous section; it usually takes place at an interaction point (either the beam dump target or collision point for accelerators), while the latter occurs in a shielded detector farther away. Furthermore, most of the searches follow simple cut-andcount strategies (up to rare exceptions, for instance missing energy searches [21,57]); we can therefore decompose the expected number of signal events as where N prod is the number of produced dark sector states, E is a detection efficiency which contains all the details about the search channel efficiency of the experiment considered, and P sig is the probability that a dark sector state leads to a signal event in the detector (for instance a decay or a scattering). In most of this section, we will place limits on our fermion portal framework by reinterpreting existing dark sector searches focused on dark matter, dark photons or dark Higgs bosons. We neglect to a first approximation the difference in kinematics. For a relativistic long-lived particle decaying through the operators (2.2) or (2.3), the detection probability depends directly on the probability of observing a decay in the decay volume of the detector.  where we have defined Γ 2 to be the decay width of the heavier unstable state χ 2 , D the distance it travels before entering the decay volume, L the distance travelled in the decay volume and γ its boost factor. We present in Table 5 the values of these parameter for various experiments. In particular, the value of the boost factor is critical in determining the lower reach of the experiments (which corresponds to the short lifetime limit for χ 2 ). We obtained the average boost factor from direct simulations of meson decays using BdNMC [58] modified to handle the decay of a dark photon mediator into two dark sector states χ 1 and χ 2 of different masses [59,60]. 11 The boost factor in the case of parton-level production (relevant for dark sector masses around the GeV) is then obtained by rescaling the average energy of the dark sector pair according to its invariant mass M 2 12 = (p 1 +p 2 ) 2 as E 2 mes − M 2 12 , with E mes the average meson kinetic energy. An important difference with respect to the fermion portal's phenomenology is that the typical fermion portal boost factor from parton level events is significantly lower than the one typically obtained from dark bremsstrahlung of dark photons (since in the latter case the dark photon carries off most of the beam energy).
Existing limits in the literature usually focus either on the case of a decaying dark photon or on the inelastic dark matter scenario, where the decay of the heaviest state is mediated by an off-shell dark photon. To estimate the sensitivity to the fermion portal we will recast these results for our effective theory in three steps: 1. We simulate the typical number of produced dark sector particles, both for the existing inelastic dark matter limits, N DP prod , (typically for δ χ ≡ |M 2 |/|M 1 | − 1 = 0.1) and for our effective theory with the required M 1 and M 2 , N eff prod . Note that we use a kinetic mixing parameter ε = 0.001 and dark sector coupling α D = 0.1 for the former, 12 and Λ/ √ g = 1 TeV for the latter (with the effective coupling's ratio either electromagnetically-aligned or Z-aligned). This choice has no consequences on the results since the scaling with respect to these parameters is trivial.
2. Focusing on the very long-lived case, where cτ γ D, L, the parameters of the experiments can be cancelled out of the ratio where g e is the effective coupling to electrons.
3. Finally, we assume that the efficiencies of the experiment between the inelastic dark matter and our effective theory case will be similar for equivalent invariant masses of the χ 1 χ 2 pair, so that using Eq. (3.1) for both the inelastic dark matter and the fermion portal case leads to This procedure can also be used to obtain the bounds for different splittings between the two hidden sector states. In this case we first proceed to estimate the efficiency around the limit by inverting Eq. (3.1). Assuming that the detection efficiency E depends dominantly on the invariant mass of the χ 1 χ 2 pair M 12 , we need to match the efficiency for the production and detection of two states with splitting δ χ to the efficiency E of two states with splitting δ χ . We have where we have assumed that the first efficiency was given as a function of M 1 , while we want the recasted one as a function of M 2 (which is the most relevant mass for large splitting). Finally, a last difficulty occurs due to the lower kinematic threshold 2m e of these decay searches based on χ 2 → χ 1 e + e − . Since most of the initial limits are estimated at small splitting, we cannot estimate the efficiency in the range M 12 ∈ [2m e , 2+δχ δχ 2m e ]. We sidestep the issue by scaling E from the available range of 12 We define the kinetic mixing as with B µν the hypercharge field strength, F µν the dark photon field strength.
where the upper limit is arbitrarily chosen to the be significantly larger than the lower threshold and nonetheless small enough so that the effect from the π 0 production threshold remains small (we have checked that the limits do not significantly depend on the precise value for the upper limit). An important subtlety, however, lies with the appearance of the lower limit in Λ, arising typically when the long-lived particle decays mainly before reaching the detector. In this case one needs to use all the geometric parameters of the experiment to estimate the decay probability. We typically use the upper limit to deduce the lower limit, using the fact that at fixed masses, the production rate scales simply as 1/Λ 4 so that the only technically challenging quantity to estimate is the decay probability ratio. More precisely, we require Assuming that the decay probability (3.2) in the lower regime is dominated by the exponentially-suppressed first contribution, the above equation leads to a simple transcendental equation on Λ low , where γ is the average boost factor of a χ 2 particle, with the relevant values collected in Table 5. These type of equations can be solved using the Lambert W -function. An expansion in logarithms in the regime of interest to us then leads to the simple expression Estimating the parameter A along with the decay probability for the upper limit P up sig implies knowing the geometric properties of the various detectors under consideration.
We have collected all the relevant parameters in Table 5. 13 Notice that this leads to conservative lower limits. Indeed, we only include the average boost factor from the dominant production mechanism, but in the short lifetime limits it is actually the high energy tail of the distribution which dominates the population in the detector and fixes the limit. Full simulation of the production and decay of the heavy states is therefore likely to lead to significantly improved lower limits with more parameter space coverage.
Most experimental limits are based on observing the creation of a positron-electron pair in the detector. The 3-body decay χ 2 → χ 1 e + e − width can be straightforwardly 13 Note that when the production of dark sector particles occurs at the LHC, the cross-section has a non-trivial dependence on the scale Λ due to the cuts described in Sec. 18 estimated as where the function G depends on the type of effective coupling to leptons. For a vector effective coupling to leptons, in the saturation limit (M 2 M 1 ) and limit of neardegenerate small splitting, we have For an axial-vector coupling, we have instead where in both cases we neglected the lepton mass in the saturation limit. In both expressions, we have assumed both M 1 , M 2 > 0. Note that as for the meson decay amplitude, in the case of a pseudo-Dirac setup (with aχ 2 γ µ γ 5 χ 1 dark fermions operator and negative mass M 1 ) one can directly use the vector coupling result with positive mass. For larger mass splitting between the dark sector states more decay channels become available. While this does not modify the expected numbers of e + e − pairs in the long lifetime limits, it can alter significantly the lower bounds from the short-lived limits. The decay into a pair of muon-antimuon opens up once ∆ χ > 2m µ , while the possible decay into hadronic states depends on the type of effective operator with quarks. The relevant channels are • Vector coupling -The dominant hadronic processes are χ 2 → χ 1 π + π − , χ 2 → χ 1 ρ and χ 2 → χ 1 ω. The corresponding decay width can be estimated with the same techniques developed in the previous section and in Appendix A, for instance with the vector meson dominance (VMD) formalism. We show the corresponding decay width in Fig. 4a.
In most cases the leptonic channels contribute significantly to the decay width so that neglecting the hadronic decay channels will be a good approximation. 14 The notable exception is the two-body χ 2 → χ 1 π 0 decay in the case of an axial-vector current. Altogether we will use the total decay width when estimating the lower limit as in Eq. 3.9. Finally, in the case of a very long-lived particle (or a completely stable one), it is possible to search for a scattering signature in the detector. However, the strong suppression from the off-shell nature of the portal implies that such limits are hardly relevant compared to the mono-X searches presented in the next section. We based our limits on the standard dark photon portal searches, such as for instance the one from Mini-BooNE [27]. We included projections for upcoming experiments based on Refs. [44,58]. 15 Since the processes involved are very similar, we use the same techniques as presented above for the very long-lived regime to recast the existing limits. Note that the coupling dependence of the limit depends on the nature of the targets and on the precise form factor for scattering off a proton or a neutron. Altogether, we approximate the scaling as being ∼ |g d | + |g u | since this gives already a qualitative understanding of the typical size of the bounds. 16

Collider and mono-X searches
Constraints from mono-photon searches at e + e − colliders are standard bounds on most models of dark photons. For the fermion portal, the off-shell nature of the mediator degrades the signal quality since one can no longer search directly for a bump in the data. A thorough analysis of this case was conducted in Ref. [65], where a limit corresponding to was found, valid when the dark sector mass is below a few GeV (at higher masses, the limited energy of the BaBar beam starts impacting the production cross-section [65]). We note that the actual limit could be enhanced with a dedicated analysis, since it was derived by Ref. [65] in a signal-only setup and with a cut-and-count approach for each bin of the reconstructed missing m χχ mass in the original analysis [66]. This implies that the most recent BaBar limits [21] using instead the final 53 fb −1 dataset do not improve the bounds straightforwardly. For Belle II, the projected limits from Ref. [65] is assuming that the dominant radiative Bhabha scattering background can be reconstructed and subtracted with a systematic uncertainty of 5%. Once again, a proper analysis including background simulation and fitting of the distribution is likely to give a stronger limit. Note that a recent study has considered displaced vertices signatures at Belle-II [67]; we leave its recasting for future work. Limits from LEP on extra-dimensions obtained at DELPHI [68] have been shown in Ref. [69] to lead to the limit We will use directly this limit, and follow Ref. [33] in adding a lower limit Λ 200 GeV to account for the breakdown of the effective theory around the LEP centre-of-mass energy.
Let us now turn to the limits at the LHC. The strongest bounds on an invisible dark sector comes from mono-jet searches, and in particular the analysis from ATLAS [70] with 36.1 fb −1 . This analysis was recasted in Refs. [33,34] for a variety of effective operators in the context of dark matter searches at LHC. We used the upper limits from these references for a vector-vector operator to extract the limit cross-section at g u = 2/3, g d = −1/3: σ lim ∼ 0.28 pb. We can then find the upper and lower limit for any effective coupling by solving (3. 16) In practice, we solve once and for all the limits in Λ as a function of 2g 2 u + g 2 d , then substituted the exact value depending on the precise values of the coefficient. This typically leads to limits of Λ > 1.3 TeV or Λ < 0.3 TeV for 2g 2 u + g 2 d = 2, with the upper limits reaching up to the tens of TeV for couplings near the perturbativity limit. Note that when considering associated production, the final limits typically turn out to grossly underestimate the reach compared to simplified models; moreover, additional effects such as multi-jet production also become relevant [71]. Our limits may then be considered a conservative estimate. 21

Invisible meson decays
Let us first discuss the limits from invisible decay of π 0 meson. The current best bound is fixed by the NA62 collaboration (see e.g. Refs. [72,73]) to be BR π 0 →νν 4.4 · 10 −9 . (3.17) Using the expression for π 0 → XX for the axial-vector effective operator one can readily deduce the limits for any given parameters. As we will see in the next sections, this limit can be quite stringent and leads to limits in the hundred of GeV on Λ.
For flavour-diagonal couplings to second and third generation quarks, the invisible decay of heavy vector mesons is constrained from BaBar and BES measurements [74,75]: More sensitive limits can be obtained from searches for invisible heavy meson decays with off-diagonal couplings. Belle has placed the strongest bounds on neutral B meson invisible decays [76]: For charged B meson invisible decays the best limits are from BaBar [77,78], though we also include a projected future bound from Belle-II [79], (3.23) BR B ± →π ± νν < 1.0 × 10 −4 (BaBar) . (3.24) Finally, we also include the following Kaon invisible decay bounds from the NA62 [73], E949 and E787 experiments [80], as well as a future projection NA62 [81].

25)
BR K + →π + a < 0.73 × 10 −10 (E949+E787) , (3.26) (BR K + →π + a < 0.01 × 10 −10 ) (NA62 projected) . (3.27) Note that for the charged Kaon decays only an invisible axion a as a final state was considered, not νν as in the other cases above, though we expect the constraints to be of similar order of magnitude. Similar limits can also be obtained for the axial-vector operator case from fully invisible decays of B and K meson (see e.g. Ref. [82]).

Astrophysical and Cosmological limits on the fermion portal 4.1 Limits from Supernova 1987A and stellar cooling
The cooling of stars and supernova 1987A are known to place strong bounds on light new physics, and as such there has been great effort to quantify these constraints (for some examples of these studies, see e.g. Refs. [18,[83][84][85][86][87][88][89][90][91][92][93][94][95][96][97][98][99][100][101]). If new light particles exist and can be produced in a star or supernova, there exist two typical regimes bracketing the bounds. In one regime the interaction strength with the SM becomes too weak so that the production of light states is no longer an effective cooling mechanism. In the other regime, the interaction with the SM becomes so strong that the light particles are produced in abundance, but are unable to exit the star/supernova. In between these two regimes is when cooling takes place (too) efficiently and can thus be excluded.

Supernova 1987A
Light particles can be produced in the proto-neutron star at the heart of a supernova. This will occur as long as the masses of the light particles are below the characteristic energy scale of the star. The core temperature of the supernova is T c ∼ 30 MeV. This enables the placing of constraints on particles with masses as large as m χ ∼ O(100) MeV, since the observed cooling of the supernova agrees within uncertainties with SM estimates.
In the case where couplings to quarks and leptons are electromagnetically aligned, assuming a thermal distribution of state near the core of the proto-neutron star, the dominant production mechanism of new light states in the supernova is through electronpositron annihilation [94]. 17 Many of the analyses in the literature consider the special case where dark fermions are coupled to a dark photon with a similar mass. Their exclusion results are therefore not straightforward to recast in the language of effective operators. However, we obtain a limit from the analysis in Ref. [94] which was performed under the assumption that the dark photon was decoupled and therefore allows for a simple recasting. The limit is obtained on the invariant mass of the pair of light fermions. Therefore, when we consider a large splitting M 2 M 1 , the upper limit on M 2 tends to be greater than naively expected. This can be seen in e.g. Fig. 5, where the splitting is large, and therefore the invariant mass of the light pair is almost entirely dominated by M 2 .
When two light states of differing mass are produced, their mass difference can have an important effect on the lower boundary of the limit on Λ/ √ g, where particles become trapped instead of streaming out of the star. For example, if the mass splitting is much greater than ∼ 30 MeV, the likelihood for the lighter χ 1 to up-scatter into χ 2 by interacting with SM particles is exponentially suppressed. If the χ 2 decay rate into χ 1 + SM is sufficiently large, this means that to a good approximation, there is no appreciable χ 2 population in the star, and there can be no annihilation or scattering of χ 1,2 to result in a trapping limit [18]. Thus in this situation, there would be no lower limit on Λ/ √ g.
If on the other hand, the decay rate of χ 2 is very small, even though the mass splitting may be too large for up-scattering to occur, there is still a significant population of χ 2 . The dark sector particles may therefore annihilate back into SM fermions, resulting in trapping and a lower limit on Λ/ √ g. The precise location of this limit would depend on the mass splitting and the decay length of χ 2 . For this reason, in our figures we show the upper limit on Λ/ √ g as a solid line, while the lower limit is dotted.
Finally, the case of an axial-vector operator is more intricate as the production is strongly modified with respect to the dark photon case. As a conservative upper limit, we thus simply use the bound from invisible π 0 decay from SN1987A cooling inferred in Ref. [102] for a core temperature of 50 MeV: BR π 0 →νν 1 · 10 −13 . (4.1) The upper limits on the effective operator are then obtained in the same way as in Sec. 3.3.
We treat the lower limit on the suppression scale in the same way as we do for the case of the vector operator above.

Stellar cooling
The characteristic temperature at the cores of Horizontal Branch and Red Giant stars is T ∼ 10 keV. This core temperature results in constraints that only apply to dark sector particles with masses as large as m χ ∼ O(50) keV. The two classes of stars differ in their densities, chemical potentials and photon plasma masses. This leads to slightly different limits being obtained from the two types of stars (see e.g. Ref. [84]). However, due to their core temperature being significantly lower than that of a supernova, the limit that can be derived on Λ/ √ g is also less strong. Indeed, in Ref. [87], it was found that the upper limit on the suppression scale of the fermion portal operator would be Λ/ √ g ∼ v, where v is the usual electroweak VEV. This can be understood intuitively from the fact that at stellar core temperatures, plasmon decay is the dominant process not only for dark sector particle production [87], but also for neutrino production. In the limit where both neutrinos and dark sector particles are massless, requiring that the luminosity in dark states not exceed the luminosity in neutrinos is equivalent to requiring that the suppression scales in their decay rates be similar in size. We do not show these bounds in our figures, since as discussed in the following section, below m χ ∼ 5 MeV, much stronger constraints can be obtained from considerations of the early universe.

Early universe and relic density bounds
Four-fermion operators like the ones considered have previously been used to describe dark matter interacting with the Standard Model in a model-independent way [33,34,65,69]. This approach for dark matter typically faces several difficulties. The final relic density obtained from the freeze-out of dark matter annihilating through a fermion portal operator typically scales as where we have considered only one vector operator of the form 1 Λ 2 (χγ µ χ)(ēγ µ e). This illustrates that the effective interactions are typically too suppressed to lead to the proper relic density through the standard freeze-out mechanism. (Note that including more effective operators and treating properly the annihilation into mesons improves somewhat the picture, as can be seen in, e.g. Ref. [35] for the case of scalar dark matter, but there is still an overabundance of thermal dark matter when Λ 200 GeV.) One can still obtain the correct relic density through the freeze-in mechanism, as pointed out in Ref. [33], although this typically implies an extremely high effective scale and adds a dependence on the reheating temperature. More generally there have been many studies of additional dynamics in the dark sector beyond the fermion portal operator which may contribute to fixing the proper relic density which require either other operators or more dark sector particles. Some of the earlier examples of such setups include Secluded DM [103] and related scenarios such as co-decaying DM (see e.g. Refs. [104,105]), or additionally Cannibal DM [106] or co-scattering [107,108], along with many more recent examples. Additionally, exotic cosmological histories can also modify the relic density, for instance through a late phase transition (see e.g. [109,110] and the subsequent literature).
The strongest limit on a possible dark matter candidate below ∼ 10 GeV comes from the cosmic microwave background (CMB) constraints on the dark matter annihilation cross-section. As shown, in, e.g. Ref. [111] unsuppressed s-wave annihilation in this case is excluded by the CMB spectral shape [112] (see also [35] for a recent study of the dim-6 case for scalar dark matter).
Finally, while we do not focus explicitly on this case in this work, very strong additional limits from the CMB arise when one of the dark sector fermion is lighter than around 5 MeV (see e.g. [113] for an up-to-date estimate). For dark fermions light enough to behave as radiation at neutrino decoupling, the strongest limits come from the effective number of relativistic degrees of freedom N eff in the early universe. As was studied in detail in Refs. [14,114,115] CMB-S4 observatories can in principle completely exclude any light relativistic relic up to arbitrarily high decoupling temperature, provided it was in thermal equilibrium with the Standard Model and that the reheating temperature was higher than its decoupling temperature. We will adapt the calculations of Refs. [14,115] to the fermion portal case in order to obtain order-of-magnitude limits. We focus on the case of a fermion portal involving electrons, 1 Λ 2 (χγ µ χ)(ēγ µ e), so that all fields involved can be considered massless in the following.
In the limit where the light particles decouple instantaneously from the thermal bath at a temperature T 0 , it will not receive the entropy released subsequently by annihilating species, so that the final effective degrees of freedom after neutrino decoupling is given by [114] ∆g * = g LS 3.38 where g LS is the number of degrees of freedom of the light relativistic relic. Translating into an effective number of neutrinos and assuming T 0 > T EW we find the lower bounds of Ref. [14] on the effective number of neutrino, In particular notice that we have kept g * (T = T 0 ) to emphasize that the bound derived in Ref. [14] is only when there are no additional degrees of freedom in the theory beyond the SM ones. This assumption does not hold by definition in our effective theory approach since we do expect new physics to occur around the scale Λ. 18 This implies that CMB-S4 experiments will not necessarily rule out every thermally coupled relic, especially in the case of a particularly rich UV sector. Current limits from the Planck experiment [112] typically also exclude light relics decoupling below the QCD phase transition at around 100 MeV. Following Ref. [14] we determine the decoupling temperature T by simply comparing the production rate of the dark sector particle χ 2 through the fermion portal operator, Γ 2 (T ), with the Hubble rate, H(T ): where we used the reduced Planck mass M pl = 2.4 · 10 18 GeV and g * (T R ) is the effective number of relativistic species at the reheating temperature T . We obtain the production rate as where we have used the equilibrium density for a Weyl fermion n eq χ = ζ(3)

4
T 3 /π 2 and introduced the thermal distribution functions f 1 , f 2 as well as a simplified Bose enhancement/Pauli blocking term P(p 1 , p 2 ) which, in our four-fermions interaction case are: (4.7) The main difference with respect to Ref. [14] comes from the dimension-6 nature of the fermion portal operator, which implies that the centre-of-mass production cross-section is of order In particular, it is a linear function of the squared centre-of-mass energy s (compared to a constant in the dimension-5 case in Ref. [14]). We then solve using the standard techniques of Ref. [116] by changing variables from p 1 , p 2 to s, E 1 +E 2 , E 1 −E 2 , including the factor P and solving numerically after extracting all T dependence, to obtain the following order-of-magnitude limits: , (4.9) 18 As an example, the full MSSM has g MSSM * = 228.75, for which the lowest value for ∆N eff is actually closer to ∆N eff > 0.01.
where the second limit depends on the reheating temperature T R since it implies that the light relics were never produced in the first place, while the first limit simply requires it to decouple prior to the QCD phase transition. Notice that the CMB-S4 limit scales as T 3/4 R as compared with T 1/2 R for the dimension-5 case [14].

Summary plots and numerical results
We present in this section some illustrative results based on the above formalism. Since the result depends strongly on the choice of effective operator, we will typically choose particular ratios g u : g d : g l and present the limits in term Λ/ √ g, where Λ is a mass scale and g represents the overall coefficient. We shall refer to g as an effective coupling though in general it will be a combination of model-dependent factors obtained by a matching calculation to a UV model. For a tree-level UV completion with O(1) couplings the scale corresponds roughly to the mediator mass. A weaker or stronger coupling could lower or raise respectively this mass scale. The usual caveats then apply regarding the regime of validity of the effective theory [117].

Vector operator
The first, straightforward case that we consider is the electromagnetically-aligned scenario which is typically obtained from an integrated-out dark photon kinetically mixed with the SM photon. The production mechanisms in this case follow closely the ones studied extensively in the literature during the last decade. Due to the off-shell nature of the process, the dark sector production at low mass is dominated by η → γχ 2 χ 1 (when it is kinematically available) and at high masses by the parton-level production. We present in Fig. 5 a summary of the current limits on this scenario as well as some projections for upcoming experiments. Limits from mono-photon signatures at BaBar and LEP, as well as the projection from Belle-II are derived following Sec. 3.2. We note that the lower limit from LEP bounds arises due to a breakdown of the effective approach around the LEP centre-of-mass energy; complete models of the dark sector including direct mediator production at LEP would then likely be excluded in this region.
The limits for LSND, CHARM, FASER, SeaQuest and SHiP that are based on the decay of long-lived dark sector states are presented for the saturation case where M 1 M 2 . The upper bound for each of those experiments can thus be seen as the maximal attainable reach. We have included current limits from LSND recasted from Ref. [60] (equivalent to the one from Ref. [118]) and the limits of CHARM by Ref. [119]. We show two future experiments as long-term prospects: a naive 10-events projection for SHiP based on the production and detection processes described above (hence not including geometric and detector efficiencies), and projected limits for FASER phase-2 (based on the study of Ref. [63]).
In all cases, the limits are recasted from a small splitting limit between M 2 and M 1 to the saturation case M 1 M 2 following the procedure presented in Sec. 3.1. The lower M 1 on the x-axis. Grey region indicates coverage from mono-photon at BaBar [65], dashed grey line Belle-II projections [65]. The shaded blue region is the mono-γ limit from LEP [69]. Limits from χ 2 → χ 1 e + e − : the green (dark green) regions are the exclusion from LSND [60] (CHARM [44]). We show a projection for FASER [63] in dashed purple and SeaQuest in red (Phase-II defined in [32]). The 10 events reach by SHiP is shown in dashed orange. The purple region represents the limit from cooling of SN1987A [93]. The normalised splitting δ χ is defined in Eq. 2.4.
limits for all these experiments are obtained following Sec. 3.1 and therefore combine two main effects: the opening of different decay channels (in particular hadronic) for larger M 2 masses, and the average boost factor associated with the various production mechanisms. As discussed before, the presented lower bounds are thus conservative since we can expect that the high-momentum tail of the dark sector particles' distribution should dominate the signal events in this region. Notice also that we did not include secondary production through up-scattering as presented in Ref. [120]. The purple region represents the limits obtained from the cooling of SN1987A, recasted from Ref. [93] as described in Sec. 4. Note that the lower limit is dashed to represent the significant uncertainty on the trapping regime.
Altogether, the existing set of limits on our fermion portal scenario presents an interesting complementarity, similar to dark photon searches, for example (but with some differences in the phenomenology, as previously discussed). Mono-X and missing energy limits tend to exclude an effective scale independently of the dark sector particles masses, but do not extend beyond Λ around a few hundred GeV. Interestingly, and contrary to the situation for a dark photon, the limit from SN1987A directly overlaps with the mono-X limits and extends to several TeV for dark sector masses below ∼ 100 MeV. Finally, the parameter space coverage of experiments based on decay searches typically extends Grey areas indicate already covered parameter space. We show prospects for FASER [63] in dashed purple and SeaQuest in red (Phase-II defined in [32]), MATHUSLA [63] in blue and SHiP in green (based on our 10 events projections).
diagonally, as could be expected from Eq. (3.10) since theses searches are most effective when the long-lived state decays a fixed distance of tens to hundreds of meters from the beam dump.
The limits we have shown in Fig. 5, for electromagnetically-aligned couplings with δ χ = 10, emphasise current limits with some representative future projections also displayed (others are omitted for clarity). In the next decade, many new experiments searching for decays of dark sector states have been proposed or are already planned. We show in Fig. 6 the projected reach for a selection of future experiments including FASER [29,63], SeaQuest (following the Phase-II proposal defined in Ref. [32]), MATHUSLA [30,63] and SHiP [31] (based on our 10 events projections), for electromagnetically-aligned couplings with splitting δ χ = 0.2.
The above results can vary depending on the effective coupling g or the mass splitting between the two dark sector states. We illustrate this dependence in Fig. 7 for the limits on Λ as a function of the effective coupling g; note that the scaling is not necessarily trivial since mono-photon and mono-jet limits from LEP and ATLAS enter above a certain energy threshold, as described in Sec. 3.2. In Fig. 8 we show how the typical limits depend strongly on the splitting in χ 2 → χ 1 e + e − decays for the CHARM and FASER phase-2 experiments (based on the limits from Ref. [63]).
As a last vector operator example for the light flavours, we present in Fig. 9  M 1 on the x-axis. Exclusion from LSND [60] (CHARM [44]) are shown in green, as well as 10 events reach by SHiP in dashed orange. The shaded blue region is the mono-γ limit from LEP [69]. The grey region is the exclusion from ATLAS mono-jet [70]. The normalised splitting δ χ is defined in Eq. 2.4. : Limits and projected sensitivity to Λ/ √ g as function of δ χ as defined in Eq. 2.4.The shaded blue region is the mono-γ limit from LEP [69]. Exclusion from CHARM [44] is shown in green, and projection from FASER [63] in purple. M 1 on the x-axis. Grey region indicates coverage from mono-photon at BaBar [65], dashed grey line Belle-II projections [65]. The shaded blue region is the mono-γ limit from LEP [69]. Limits from χ 2 → χ 1 e + e − : the dark green region is the exclusion from CHARM [44]. We show a projection for FASER [63] in dashed purple and SeaQuest in red (Phase-II defined in [32]). The 10 events reach by SHiP is shown in dashed orange. The purple region represents the limit from cooling of SN1987A [93]. The normalised splitting δ χ is defined in Eq. 2.4.
in the case of proto-phobic couplings, The main difference with the previous electromagnetically-aligned case is that π 0 decay production is strongly suppressed, so that searches at the LSND experiment do not lead to a significant limit. On the other hand, other beam dump experiments rely mostly on the η meson decay which is only mildly modified, as can be seen in Table 2. Fig. 9 is for δ χ = 10. While most of the limits above are based on the coupling with first generation fermions, effective vector operators for heavy flavours can also be constrained using the same techniques. We present the limits for SHiP based on K and B meson decays in Fig. 10. The invisible decay bounds are also shown as labelled for BaBar, Belle (II), NA62 and E949/787. We see that the heavier masses of the mesons involved can significantly extend the fermion portal sensitivity to higher effective scales. Notice that the invisible meson decay constraints appear to exclude a large region of the parameter space that SHiP will probe, in particular for the K meson production case, but we recall that we assumed here the same scale suppression and couplings for both production and decay.

Axial-vector operator
As a first example of limits based on the axial-vector operator, we focus on the Z-aligned limit with We present the result of current limits and a representative selection of future sensitivities, for a splitting δ χ = 10 in Fig. 11. Notice that the two-body meson decay production mechanism strongly enhances the limits at low masses. In particular, the recasting of LSND searches leads to bounds up to the TeV scale in this case. An interesting feature of the lower limits for χ 2 decay is the strong reduction of the limits for M 2 − M 1 > M π due to the opening up of the χ 2 → χ 1 π 0 decay channel, as described in Sec. 3.1. In Fig. 12 we emphasise the projected reach for future experiments, including FASER [63], SeaQuest (following the Phase-II proposal defined in [32]), MATHUSLA [63] and SHiP (based on our 10 events projections), again for the Z-aligned couplings but with a smaller splitting δ χ = 0.2. Shown also shaded in grey are the present exclusions from SN1987A [102], LSND [60], CHARM [44], LEP [69], BaBar [65] and NA62 [73].
Finally, in Fig. 13, we consider the case of an axial-vector effective operator with "universal" couplings g u : g d : g l = 1 : 1 : 1 . M 1 on the x-axis. Grey flat regions indicates mono-γ limit at BaBar [65], dashed grey line Belle-II projections [65]. The light grey regions at low M 2 indicate NA62 π 0 →inv limits [73]. The shaded blue region is the mono-γ limit from LEP [69]. Limits from χ 2 → χ 1 e + e − : the green (dark green) regions are the exclusion from LSND [60] (CHARM [44]). We show a projection for FASER [63] (dashed purple) and SeaQuest in red (Phase-II defined in [32]). The reach of SHiP for 10 signal events is shown in dashed orange. The purple region represents the π 0 → νν limit from [102] based on cooling of SN1987A. The normalised splitting δ χ is defined in Eq. 2.4.  Figure 12: Prospective sensitivity to the axial-vector operator effective scale Λ/ √ g from future intensity experiments in the case of Z-aligned effective couplings as a function of M 2 on the x-axis for δ χ = 0.2. Grey areas indicate already covered parameter space. We show prospects for FASER [63] in dashed purple and SeaQuest in red (Phase-II defined in [32]), MATHUSLA [63] in blue and SHiP in green (based on our 10 events projections). The normalised splitting δ χ is defined in Eq. 2.4.  [65], dashed grey line Belle-II projections [65]. The light grey regions at low M 2 indicate the exclusion from NA62 π 0 →inv limits [73]. The shaded blue region is the mono-γ limit from LEP [69]. Limits from χ 2 → χ 1 e + e − : the dark green region is the exclusion from CHARM [44]. We show a projection for FASER [63] in dashed purple and SeaQuest in red (Phase-II defined in [32]). The 10 events reach by SHiP is shown in dashed orange. The SN1987A limit is not shown as explained in the text. The normalised splitting δ χ is defined in Eq. 2.4.
This translates into an effective pion-phobic scenario due to the form of the effective coupling to pions, as shown in Table 2. As in the proton-phobic case of the last section, the limit from LSND vanishes. Additionally the heavy state decay χ 2 to a π 0 χ 1 is strongly suppressed, extending downwards the limits. Finally, the limits from SN1987A obtained by considering the invisible decay of neutral pions [102] no longer apply. A limit from SN1987A from e + e − annihilation in the supernova core should still apply, but we have not computed it here. It would likely be similar to the e + e − annihilation constraint for the vector operator shown in e.g. Fig. 5.

Small mass splitting
In most of the plots presented above, the limits coming from the decay of the heavier χ 2 dark sector state were important but highly dependent on the lifetime. Since the lifetime of the heavy state scales as ∆ 5 χ in the limit of small splitting between the dark sector masses, these limits are reduced for small splitting. We illustrate this aspect in Fig. 14 on the top and bottom plots for the (Z-aligned) axial-vector and (electromagneticallyaligned) vector case with δ χ = 0.01 and 0.05, respectively. In particular, we have rep- Grey flat regions indicates coverage from mono-photon at BaBar [65], dashed grey line Belle-II projections [65]. The light grey regions at low M 2 indicate the exclusion from NA62 π 0 →inv limits [73]. The shaded blue region is the mono-γ limit from LEP [69]. Limits from χ 2 → χ 1 e + e − : the dark green region is the exclusion from CHARM [44]. The 10 events reach by SHiP is shown in dashed orange. Scattering limits at MiniBooNE [27] are shown in light green, and SBND [58] and NOνA [44] in orange and red. The normalised splitting δ χ is defined in Eq. 2.4. For limits from χ 2 → χ 1 e + e − , the dark green region is the exclusion from CHARM [44] and the 10 events reach by SHiP is shown in dashed orange. We also show a projection for FASER from [63]. The normalised splitting δ χ is defined in Eq. 2.4.
resented the limits based on scattering of the lighter state in MiniBooNE (based on Ref. [27]), SBND (projections from Ref. [58]), and NOνA (from the 3 · 10 20 protons on target projection of Ref. [44]). While these limits are not competitive with the missing energy searches from BaBar, they may become relevant in a more lepton-phobic scenario. Notice additionally that the lower mass thresholds for the χ 2 → χ 1 e + e − decay to be allowed are shifted to higher masses relative to figures where the mass splitting is large.

Concrete scenario: GeV scale dark photon
We end with a practical application of the approach presented above, in terms of the familiar dark photon benchmark as a possible UV completion of the fermion portal. Limits on dark photons decaying invisibly to dark sector fermions with a couplings g D are currently relatively weak, in the tens of GeV range, with the current best limits arising from LEP [121]. Interestingly, such a heavy dark photon V has a sizeable mixing with the Standard Model Z boson, leading naturally to an axial vector current coupling with the SM fermions. One can go from the effective approach to the simplied model (in the off-shell dark photon limit) using Noting that at first order in the kinetic mixing ε the axial vector coupling are δ 2 - where c W is the cosine of the Weinberg angle. The vector couplings are not significantly modified as long as δ 2 1: An important point is that while production of dark sector states can now proceed through either of the operators (including the significant boost observed at low masses in the axial-vector case), the lifetime depends on the sum of the two decay width. In particular, "mixed" contributions where production proceeds through the axial-vector operator and decays through the vector one dominates at low masses (when meson decay production for π 0 and η dominate, as seen in Sec. 2). We illustrate these effects and summarise the bounds in Fig. 15, where we show the limits on ε for M V = 20 GeV. Note that the limit from FASER are conservative in that this experiment will have access to enough energy to produce on-shell a dark photon of such mass.

Conclusion
Light dark sectors are a class of new physics beyond the SM that present special challenges and opportunities for discovery at the intensity frontier. Since they are neutral under the SM gauge groups, dark sector fields only interact with the SM through so-called portal operators -gauge singlet combinations of SM fields. Much work has been done on studying the phenomenology of the three lowest-dimensional portal operators: the Higgs, vector, and neutrino portals. Here we have presented a study of the next lowestdimensional "fermion" portal. We focused on the dimension-6 four-fermion operator combination of a fermion portal to a pair of light dark sector fermions. The higherdimensional nature of this portal can arise naturally from the off-shell limit of one of the renormalisable portal interactions. The scale suppression could moreover explain the weakness of the interaction between visible and dark sectors. Our effective field theory approach has the advantage of encapsulating the phenomenology of light dark sectors interacting with the Standard Model through heavy mediators in full generality.
In our study the typical production mechanisms at intensity frontier experiments present several interesting modifications compared to the standard vector or scalar portals, for instance. Typically, light meson production is subdominant due to the scale suppression arising from the higher-dimensional nature of the portal. An interesting exception is that light pseudo-scalar mesons can have a fully invisible two-body decay in the case of the axial-vector fermion portal, leading to a strong enhancement of the reach of low-energy experiments for such scenarios. Additionally, different phenomenology may be exhibited in bremsstrahlung production of light fermion pairs via an off-shell mediator at electron beam dump experiments. Despite the scale suppression, this may be a relevant mechanism that we will consider in future work.
Similarly, the detection strategies adopted for the renormalisable portals are modified by the off-shell nature of the fermion portal. In particular, missing energy limits (e.g. mono-photon searches at BaBar or Belle-II) only lead to weak bounds in our case due to the absence of a resonance in the missing mass spectrum. Limits from the scattering or decay of heavy dark sector states can also be adapted to the case of a fermion portal. To place most of our limits on the fermion portal, we simulate the production and decay of light states and use the comprehensive existing analyses for renormalisable portals, such as dark photon and inelastic dark matter models, to estimate experimental efficiencies. Note that in this analysis, we do not consider possible differences in experimental efficiencies due to kinematics. We also do not account for possible renormalisation group running and consequent mixing of fermion portal operators between the various scales involved, which can vary from MeV for decay processes to TeV for parton-level production at the LHC (see e.g. Refs. [122,123] for the dark matter case). Our framework for applying limits on the fermion portal is available as the public DarkEFT code (described in Appendix B).
We outlined the parameter space coverage of existing and future experiments for both vector and axial-vector operators, considering flavour scenarios such as electromagneticallyaligned or Z-aligned couplings and proto-phobic or pion-phobic couplings. Depending on the nature of the operator, we found that the effective scale for light flavours was typically constrained to lie in the hundreds of GeV to TeV range for effective couplings g ∼ O(1). The combination of missing energy searches typically lead to the requirement that Λ/ √ g 500 GeV for dark fermion masses of up to a few GeV. There is a gap in the limits where the effective field theory approach breaks down when considering searches at LEP. This would typically be excluded in a complete model including a kinematically accessible mediator. This gap can also be mostly covered from limits from invisible meson decays when those are available. For smaller dark sector masses up to m χ ∼ 100 MeV, astrophysical constraints from SN1987A cooling play a key role in constraining the effective scale up to several TeV, although the lower "trapping" bound has a strong model-dependence as discussed in Section 4. Limits from experiments involving heavy B and K mesons could extend the sensitivity to effective scales in the tens of TeV. While we focused mostly on dark sector fields heavier than several MeV, we also derived an order of magnitude estimate of around 5 TeV from Planck on the limits on the effective scale from N eff when the dark sector fields behave as relativistic matter. This paper has focused on fermion portal operators which are either flavour-diagonal in both lepton and quark sectors or flavour-breaking in the quark sector. In particular we do not currently consider flavour-breaking operators in the leptonic sectors and leave to future work limits on the neutrino-based operators. These should typically be generated in the UV along with leptonic ones -especially when considering the axial-vector fermion portal operator. We also note that the restriction to fermion portal operators to a pair of dark sector fermions, with a (χΓχ) structure in the dark sector, is actually not particularly restrictive for characterising the fermion portal more generally, since most of the phenomenology will depend on the SM part of the higher-dimensional operator. Our results can then be considered as good guidelines for fermion portals to other dark sector combinations (for instance S∂ µ S with S a dark scalar).
For future work there are more directions where the bounds can be either refined or improved, some of which were discussed above. One important addition would be to include and simulate the production rates from heavy meson decay and their detection prospects in terms of decay or scattering in high energy frontier experiments. The order of magnitude estimate for SHiP presented here points to limits significantly stronger than those from standard light invisible meson decays. Another refinement of the limits would be to simulate more thoroughly the production and decay of heavy dark sector states through the fermion portal; in particular near short lifetime limits, where our conservative estimates could be improved since the high-energy tail of the spectrum will dominate the expected events. Other portal operators could also be investigated.
In this work we have taken a step towards a systematic study of the Standard Model portal operators through which dark sectors necessarily interact. As the next decade begins, the intensity frontier of particle physics could enable a thorough exploration of the universe's dark sectors.
for the case of a vector meson V with polarisation vector ε λ . The amplitude then follows straightforwardly. For ρ, ω, we use directly the summary from [124], which accounts for the mixing effects between vector mesons. The decay constants associated to each quarks are: We can therefore straightforwardly project our operators on the interpolating current for the vector mesons to obtain the effective couplings presented in Table 2.
For the direct decay of a pseudo-scalar meson M → χχ in the case of axial-vector current, we follow the approach for neutrino decays presented in [41], and extend it to the η and η case using [125] to account for the mixing effects. In more details, we introduce the currents A µ 0 , A µ 8 as: And the corresponding decay constant f 8 and f 0 defined as In the two angle mixing scheme (see e.g [125]), the η and η are then represented by: where θ 8 ∼ −22 • and θ 8 ∼ −9 • . Finally, we can project our set of operators on the A µ 0 , A µ 8 basis to obtain (using the usual shorthand notation for cos and sin as c, s): . (A.9) Using the fitted values f 8 = 1.28f π , f 0 = 1.2f π from [125], we obtain the results of Table 2. Light meson three-body decays In the case of an "associated" radiative decay M → Xγ. We obtained the results from the main text using two distinct procedures. The first approach relies on the Vector Meson Dominance (VMD) advocated in [126] which has already been used extensively for the case of new dark vector particles searches (see e.g. [64,127]). We define the U (3) generators for the relevant mesons as: where we have used the same simplified approach with a single η − η mixing angle θ with cos θ ∼ √ 6/3 and sin θ ∼ −1/3 as in [64,127] (see also [128]), based on [129]. Furthermore, we can define the electromagnetic and dark coupling matrices as: Using the Feynman rules following [126] and defining the coupling g ρππ ∼ 6.1, we first consider the amplitude for the process V → χχ, where V is a vector meson of polarisation vector V : This leads to effective coupling constants within ∼ 10% of the complete results from [124] used in Table 2 (this result can also be recovered by taking the heavy dark vector limit based on the expressions from [64,127]). Based on the calculation of the π 0 → γγ amplitude from [126], we can now estimate the amplitude corresponding to the threebody decay of a pseudo-scalar meson P (p π ) → γ(q)χ 1 χ 2 from Fig. 16a as: × ε µνρσ q µ ε ν (p π σ − q σ )(ū(χ 1 )γ ρ v(χ 2 )) where we have used g ρωπ = − 3g 2 ρππ 8π 2 Fπ , with F π 93 MeV. Note that we have neglected the momentum dependence in the vector propagator so that it simply amounts to inserting a 1/M 2 V factor; we briefly discuss this approximation at the end of this Appendix. We can deduce the effective coupling g P from the first line of this equation. By comparing with the coupling in the π 0 → γγ case, we obtain: (A. 13) This expression agrees with the one of Table 2 (based on the second approach below) at 10% level for g u , g d and 30% level for g s , all the observed discrepancies are related to the single angle scheme, and replacing the T P of η and η by the two mixing angles scheme which we will use below leads to perfect agreement. The second approach builds directly on the effective amplitude for the radiative decay of a π 0 into a photon with polarisation ε λ and a dark vector boson V of polarisationε κ : A P γV ⊃ eg Q 4π 2 F π ε µνρσ q µ ε νε κ (p P σ − q σ ) , (A.14) where the effective coupling is estimated from the triangular diagram presented in Fig. 16b as (see [127,129,130]) 15) and where the T P matrices arises according to the quark content of the interpolating currents of the pseudo-scalar mesons. We obtain the amplitude in our case for a vectorcurrent interaction by replacing the dark photon polarisation vector by the dark vector current (χ 1 γ ρ χ 2 )/Λ 2 (which is of course similar to integrating out the dark photon): A P γχχ ⊃ eg Q 4π 2 F π Λ 2 ε µνρσ q µ ε ν (ū(χ 1 )γ ρ v(χ 2 ))(p P σ − q σ ) . (A. 16) While the matrices T P have been described above, for the case of η and η meson, we can go beyond the single mixing angle approach presented above using [125]. Defining

B Available limits and DarkEFT companion code
In order to simplify the use of the results presented in this paper, we have released a companion code DarkEFT, written in python and available at: https://github.com/Luc-Darme/DarkEFT. Its main features are: • A database of relevant analysis and limits, along with the relevant references and a small description.
• Amplitudes for various relevant production and decay mechanisms for dark sector states within the effective field theory presented above.
• A set of tools to recast the stored limits to the fermion portal case.
Importing the main python module of the code can done as import L i m i t s L i s t a s l i m DarkEFT allows to very simply recast existing limit for any choice of the effective couplings. As an example, recasting the limits from the MiniBooNE collaboration from light dark matter scattering presented in [27], for a electromagnetically-aligned vector operator and for a splitting of 25% between χ 1 and χ 2 can be done by g e f f e m={" gu11 " : 2 / 3 . , " gd11 " : − 1 / 3 . , " g l 1 1 " : −1. The full details and possibility of the code are presented directly in the Readme file. The current sets of implemented limits are presented in Table 7.
Notice that while most of the limits are obtained from recasting existing works, the tools presented in this article can also be used to obtain naive estimate for various setups (not including the detection and geometric efficiency). for instance, limits from long-lived state at SHiP are in fact a 10 events line obtained in this way.
The code also allows to recast a large sets of experiments and print directly the results to files. For instance the following code  Table 7: List of experimental searches currently implemented, along with some details about the process. Note that for SHiP decay limits, the 10-event limits are obtained directly using the tools described in this paper E x p e r i m e n t s L i s t=np . a r r a y ( [ " l s n d d e c a y " , " charm decay " , \ " babar monogam " , " belle2 monogam " ] ) g e f f Z a l ={" gu11 " : 1 / 2 . , " gd11 " : − 1 / 2 . , " gd22 " : − 1 / 2 . , " g l 1 1 " : −1/2 , " g l 2 2 " : −1/2} Lim , L a b e l L i m i t = l i m . GetLimits ( E x p e r i m e n t s L i s t , 1 0 , g e f f Z a l , "AV" , True ) creates a list of analysis to be recasted, loads Z-aligned couplings for an axial-vector effective operator and prints the resulting recasting into files. The output variable Lim contains a python dictionary of array (M 2 , Λ) with keys given in ExperimentsList.
We have finally enclosed with the distributed code various example files along with some plotting routines, including those for generating all the plots in this paper.