Testing freeze-in with axial and vector $Z'$ bosons

The freeze-in production of Feebly Interacting Massive Particle (FIMP) dark matter in the early universe is an appealing alternative to the well-known - and constrained - Weakly Interacting Massive Particle (WIMP) paradigm. Although challenging, the phenomenology of FIMP dark matter has been receiving growing attention and is possible in a few scenarios. In this work, we contribute to this endeavor by considering a $Z^\prime$ portal to fermionic dark matter, with the $Z^\prime$ having both vector and axial couplings and a mass ranging from MeV up to PeV. We evaluate the bounds on both freeze-in and freeze-out from direct detection, atomic parity violation, leptonic anomalous magnetic moments, neutrino-electron scattering, collider, and beam dump experiments. We show that FIMPs can already be tested by most of these experiments in a complementary way, whereas WIMPs are especially viable in the $Z^\prime$ low mass regime, in addition to the $Z^\prime$ resonance region. We also discuss the role of the axial couplings of $Z^\prime$ in our results. We therefore hope to motivate specific realizations of this model in the context of FIMPs, as well as searches for these elusive dark matter candidates.


Introduction
One of the most intriguing puzzles of Particle Physics and Cosmology is the nature of dark matter (DM), a non-relativistic matter component that makes up about 27% of the current cosmic energy [1]. The observational evidence for the existence of DM is overwhelming, relying on its gravitational interaction with ordinary matter. DM is required to explain the anisotropies of the Cosmic Microwave Background (CMB), the flatness of galaxy rotation curves, and the large-scale structure of the Universe. Nevertheless, despite the large number of DM candidates that emerge from theories beyond the Standard Model (SM), the nature of DM remains unknown (see Ref. [2] for a review). Moreover, although many experimental searches have been carried out over the past decades, DM has evaded detection.
Among the DM candidates, weakly interacting massive particles (WIMPs) are certainly the most popular ones. These particles were kept in thermal equilibrium with the SM bath in the early Universe, and just when their interaction rates fell below the Hubble rate -meaning that the interactions could not keep up with the expansion of the Universe -they decoupled from the cosmic bath and its abundance froze-out, yielding the observed DM relic density [3]. Since the cross-sections that generate the observed WIMP abundance are typically at the electroweak scale, currently being probed at direct, indirect detection and collider experiments, WIMP models deserve to be completely explored. Although WIMPs dominate the searches for dark matter and might soon be discovered, the current situation of stringent constraints [4,5] motivate us to look for alternative scenarios.
The special case of a pure vector Z which kinetically mixes to photons, the dark photon, has been extensively studied in the literature, both in the context of WIMPs [56][57][58][59][60][61] and FIMPs [14,16,19,21]. It is already known that sub-GeV dark photons have the interesting feature of enhancing direct detection rates and rendering FIMPs testable [19].
In this work, we intend to further investigate the phenomenology of FIMPs in a Z portal model as well as access the viable WIMP parameter space in a wider Z mass range, from MeV to PeV. We therefore study the freeze-out and the freeze-in production of a Dirac fermion DM candidate, χ, which only interacts with the SM through a Z boson. We consider both vector and axial-vector couplings in the Z currents, for a wide range of values. We evaluate how the region of our parameter space providing the right amount of DM is constrained by direct detection (XENON1T), colliders (LHCb, ATLAS, LEP II, BaBar), neutrino-electron scattering, atomic parity violation, electron and muon anomalous magnetic moments, and electron beam dump experiments (E137, E141).
In the presence of axial couplings, direct detection bounds are easily evaded and would hardly probe the FIMP scenario. Nevertheless, we show that most of the experiments considered in this work are already testing the parameter space of FIMPs. As a consequence, even though axial Z 's are tightly constrained, especially in their light regime [62,63], they provide a viable framework for FIMP dark matter.
The paper is organized as follows: in Section 2, we introduce our model. In Section 3, we describe how the present DM abundance can be achieved both in the context of the freezeout and the freeze-in mechanisms, whereas in Section 4 we show how the parameter space can be constrained using information from various experiments, as well as cosmological and theoretical bounds. Finally, in Section 5 we present our main results and in Section 6, we conclude.

The model
In this work, we address the possibility of testing feebly interacting dark matter whose interactions with ordinary matter are mediated by a Z , the massive gauge boson coming from an extra gauge U (1) symmetry, U (1) .
We consider an extra Dirac fermion, χ, as our dark matter candidate. When both the SM fermions (f ) and χ are charged under U (1) , the relevant Lagrangian is given by where m χ and m Z are the dark matter and Z masses, and V χ,f and A χ,f are respectively vector and axial dimensionless couplings. These are the six free parameters considered in our analysis. In our convention, in terms of the dark gauge coupling and chiral charges, we have Because SM neutrinos are left-handed, we will always fix V ν = A ν . For simplicity, we neglect possible mass or kinetic mixing between the SM hypercharge gauge boson and the Z . Note that kinetic mixing can also be generated at loop level via SM fermions. However, given that the kinetic mixing has been constrained to be small [64], these small corrections do not qualitatively change our results.
As we consider in Section 4, many experiments search for dark matter and Z 's. In order to evade the current stringent bounds while still hoping to discover them in upcoming experiments, many realizations of a Z portal were proposed.
Direct detection experiments are less sensitive to dark matter in certain Z models. For purely axial Z 's, with A f = 0 and V f = 0 (which means X f L = −X f R for all fermions) [40,42,[65][66][67], the scattering between dark matter and quarks is mainly spin-dependent, which is less constrained 1 . Furthermore, if we had chosen χ as a Majorana fermion (neutral), its vector current, associated with its charge under U (1) , would have to be exactly zero. It is also interesting to notice that potential signals of new physics can be explained by considering non-vanishing axial couplings, as in the case of MeV anomalies [68] and the Galactic Center gamma-ray excess [40].
The constraints on Z interactions with standard fermions are also quite stringent. Flavour and generation specific realizations are then usually invoked. This is the case of leptophilic models, with the Z coupling mainly to leptons [69][70][71][72][73][74][75][76][77][78]; leptophobic models, where the Z couples mainly to quarks [29,76,[79][80][81][82]; and models where the Z couples to the SM fermions in a non-universal way, for instance, with couplings to only the third generation quarks or first generation leptons [40,45]. Allowing Z 's to decay predominantly into invisible states (in our case, V χ , A χ V f , A f ) also makes collider bounds weaker. Other possibilities are sequential Z models, where V f and A f are the same as the SM Z bosons couplings, or simply re-scaled [35,43,82,83]; and kinetic mixing portals [84-  89], coming from the fact that the kinetic terms of two U (1) gauge bosons are in general non-diagonal. The case of a dark photon (A f = 0, or X f L = X f R , for all SM fermions), typically kinetically mixed to photons, has been extensively studied in the literature and is often regarded as a target for future experiments [14,[57][58][59][60][61]. To motivate model-dependent analyses with specific charge assignments in the context of FIMPs, in what follows we focus on universal V f and A f couplings. Finally, we would like to stress that realistic and UV complete Z models with a low-energy Lagrangian as in Eq. (2.1) requires the introduction of additional fields. In the presence of axial couplings, gauge invariance of the Yukawa sector requires the introduction of additional scalars [68]. On the other hand, additional scalars might also be needed to generate the masses of the χ and Z . Moreover, the presence of axial currents introduce triangle anomalies which must be cancelled. This is typically done by invoking new fermions charged under both standard and dark U (1)'s [40,66,90]. To keep our analysis as model-independent as possible, we assume that all these additional fields are too heavy to impact the production of dark matter 2 . Our results would therefore be applicable to any chiral charge assignments under these assumptions.
In the next section, we show how both the strength of the couplings and the masses of χ and Z determine the way dark matter is produced in the early universe.

Relic density
Our dark matter candidate, χ, is produced in the early universe either through the freeze-out or freeze-in mechanism, depending on whether or not it has achieved equilibrium with the thermal bath species, respectively. The processes that change the number density of DM are Z decays and annihilations through t and u-channels, and SM fermion annihilations through s-channel, as well as their backreactions in the case of freeze-out. These are illustrated in Fig. 1. In order to determine the DM relic density today, we must solve the Boltzmann fluid equation for its number density, n χ . We solve it in terms of the DM yield Y χ = n χ /s, with s the entropy density. Taking into account the leading interactions in our model, we have, where x ≡ m Z /T , T is the temperature of the thermal bath, H is the Hubble rate, g * s (x) = 1 + 1 3 d log gs(T ) d log T , with g s (T ) the number of effective degrees of freedom associated with entropy, Y eq χ is the DM yield at equilibrium, and the γ's are the reaction rate densities for each process, given below 3 .
The contribution of Z decays and annihilations to the production of χ depends strongly on whether they are part of the SM thermal bath, as otherwise they would be underabundant compared to the SM fermions.
For the decay process, we have where K n is the modified Bessel function of the second kind of order n, n eq i is the equilibrium number density of species i, and the decay rate of Z into DM χ is given by, For the self-annihilation processes of bath species b, we have γ bb→χχ (x) = (n eq b (x)) 2 σv bb→χχ = (n eq χ (x)) 2 σv χχ→bb where σv are the thermally averaged cross sections. In the expression above, the approximation holds for initial states obeying Maxwell-Boltzmann statistics, and |M| 2 is the nonaveraged squared amplitude. For the Z Z →χχ process, at high temperatures (T >> m Z ), we have γ Z Z →χχ ∝ T 6 . The s-channel rates exhibit a resonance regime, in which we can use the narrow width approximation, and off-shell regimes for light Z (m Z T, √ s) and heavy Z (m Z T, √ s). The explicit expressions for the rate densities of our 2 → 2 processes, as well as approximations, can be found in Ref. [92].  Figure 2: Ratios of the production rates (γ/n) to the Hubble rate (H) for the three DM production processes: Z annihilations (in orange) and decays (in brown), and f annihilations (in blue). The dotted black line shows where γ/n = H. For this set of parameters, we can see that the freeze-in regime is achieved for m Z = 10 3 GeV (dashed curves), whereas freeze-out occurs for m Z = 10 −3 GeV (solid curves).
Roughly speaking, if all the reaction rates are always smaller than the cosmic expansion rate (γ/n b H), χ would not be able to thermalize and freeze-in production is then possible. Otherwise, freeze-out would take place, as n χ becomes comparable to its equilibrium value and the backreactions begin to be relevant in Eq. (3.1).
In Fig. 2, we show the ratio of the reaction rates to the Hubble rate as a function of T , for the Z decays (in brown) and annihilations (in orange), and SM fermion annihilations (in blue). The dotted horizontal line, where γ/n b = H, indicates where we would have a rough boundary between freeze-in and freeze-out. Of course, the stronger the overall coupling, which is a function of the vector and axial couplings and the masses, the easier it is for χ to be produced via freeze-out. We illustrate the impact of m Z on the production regime by considering it to be in the MeV scale (solid curves) and in the TeV scale (dashed curves), for a given set of couplings and DM mass as indicated in the figure. Note that a solid brown curve is not present since the decay is not kinematically allowed. As we can see, even for tiny couplings between χ and Z (V χ = A χ = 10 −10 ), the exchange of a light enough Z is able to sufficiently enhance the s-channel cross section to thermalize DM. On the other hand, the Z does not need to be much heavier than χ to suppress the rates and enable the freeze-in regime.
For each point in our parameter space, comprised of m χ , m Z , V χ , A χ , V f , and A f (assuming universal couplings to SM fermions), we check whether the sum of all the kinematically available processes are sub-Hubble for temperatures above their Boltzmann suppression. We therefore have in the next figures dashed curves indicating the regime boundary, between freeze-in and freeze-out. Furthermore, we do the same analysis for the processes which would thermalize Z with the SM thermal bath. The region in our parameter space where Z 's are not part of the SM bath will be indicated in blue and labeled "non-thermal Z ".
In Fig. 3, we show in the plane (m Z , A f ) the contours of observed relic density of χ as inferred by Planck, Ω 0 χ h 2 0.12 [1] (solid curves), produced either by freeze-out or freezein, to the left and to the right of the regime boundary, respectively. The regions between the freeze-out and freeze-in contours would lead to an overproduction of dark matter, which would overclose the universe, and are therefore excluded by Planck. As we have said, the purely vector Z case is extensively explored in the literature, so in this figure we show how different m χ values impact the contours in the purely axial Z case. In magenta (green) we see the DM relic density contours as solid curves and the regime boundaries as dashed curves, for m χ = 1 GeV (m χ = 100 GeV). Also, we have chosen a small value for the dark matter coupling (A χ = 10 −10 ) in order to have a better understanding of the contours in the freeze-in region.
As expected, we see in Fig. 3 that the region of our parameter space in which the Z is not part of the SM thermal bath (blue region) corresponds to small couplings, in this case, A f . However, the thermalization between the Z and f 's depends strongly on m Z . The annihilationsf f → Z Z dominate the Z thermalization for lower masses, with rates increasing rapidly as the Z becomes lighter. In contrast, inverse decaysf f → Z dominate for larger masses, with rates directly proportional to m Z . In the region of a non-thermal Z , only the s-channels contribute to the DM relic abundance, since n Z n eq Z ∼ n eq f 4 . In this case, the observed DM abundance can only be achieved for m Z > 2m χ . Let us now focus on the regime boundary, indicated by dashed curves in Fig. 3 (as well as in Fig. 4). As stated above, this boundary lies where the sum of the reaction rates of the processes depicted in Fig. 1 equals the Hubble rate. Z annihilations, which are independent of the SM couplings, dominate when m Z < m χ . We therefore recognize that the vertical part of the regime boundaries is due to t/u-channels. For larger values of SM couplings, the SM fermion annihilations start dominating the thermalization of χ, corresponding to the diagonal parts of the regime boundaries. Decays and inverse decays can only dominate the χ thermalization for much heavier Z 's, not considered in the parameter space of interest in this work.
As we can see in Fig. 3, the lighter the DM, the lighter the Z must be to provide the observed DM abundance. Accordingly, the regime boundary also shifts with the DM mass. The features of the freeze-out and freeze-in contours are discussed in the following subsections.

Freeze-out regime
In the freeze-out regime, dark matter was initially part of the thermal bath, so the initial condition for Eq.
The final relic density is found by using the usual freezeout approximation [3].
In Fig. 3, the contours providing the observed abundance of DM today are seen in the left upper corner (solid curves above the regime boundaries). In this case, WIMPs annihilate predominantly into SM fermions. As we will see later, though, WIMP annihilation into Z 's dominate for smaller values of m Z . Since the case of WIMP dark matter with a Z portal is well-known, in Fig. 3 we choose to focus on the behavior of the FIMP relic density contours. Had we chosen much larger values of A χ (see for instance Ref. [42]), as we do in Section 5, the parameter space providing viable WIMPs would be larger.

Freeze-in regime
In the freeze-in regime, dark matter is assumed to be initially absent in the early universe, so the initial condition for Eq. (3.1) is Y χ = 0. Since Y χ is always much smaller than Y eq χ at least until the production finishes, we can safely neglect the term with Y 2 χ in the right-hand side of Eq. (3.1). We therefore simply integrate Eq. (3.1) from the reheat temperature T R , which we take to be 10 14 GeV for the entire analysis, down to the current temperature.
The freeze-in process can finish at the lowest scale available (infrared freeze-in), just like the usual case of freeze-out, or at the highest scale (ultraviolet freeze-in), in which case the relic abundance depends on the reheat temperature. Processes whose production rate densities have a high enough temperature-dependence can lead to ultraviolet freeze-in (as in the non-resonant Z portal of Ref. [93]). In the radiation era, this happens if the main process has γ ∝ T n with n > 5. Ultraviolet freeze-in is achieved via higher dimensional operators. In our model, we have two of such processes: the heavy Z regime of the s-channel, which happens as a four-fermion interaction and features γf f →χχ ∝ T 8 /m 4 Z , and the t/u-channel annihilation of Z at high temperature in the presence of axial couplings (see discussion in Therefore, in regions of the parameter space where these processes dominate, all FIMPs were produced around T R . The Z →χχ decay process freezes-in at T ≈ m Z , when the number density of Z becomes Boltzmannsuppressed. Thef f →χχ process also freezes-in at T ≈ m Z whenever the resonance is allowed. In the case that m Z < max[m χ , m f ], the process becomes Boltzmann-suppressed before the resonance occurs, therefore freeze-in happens at the smallest scale kinematically The contours of observed abundance produced through freeze-in shown in Fig. 3 have interesting features due to the different production processes dominating freeze-in, as described below.
Z annihilation This process is only relevant for the achievement of the correct relic density in the presence of axial couplings, in which case it leads to UV freeze-in. In the high energy limit, when the Z momentum is much higher than its mass, the contribution of Z annihilations to the relic density is given by where for simplicity we set all degrees of freedom constant, g ef f ≡ g s = g e . For m χ = 1 GeV (m χ = 100 GeV), the t/u-channel sets the correct relic density for m Z ∼ 0.2 GeV (m Z ∼ 5 GeV).
f annihilation in the light Z regime For larger values of A f , the s-channels begin to dominate, changing the slope of the contours. While m Z < 2m χ , the s-channels are in the light Z regime (s m 2 Z ) and their contribution to the relic density is found to be (3.6) As we can see from this expression, this process leads to a relic contour independent of m Z in the absence of axial couplings (see top right panel of Fig. 4 below).
f annihilation in the resonant Z regime / Z decay When m Z > 2m χ , the on-shell production of Z , and the subsequent Z decay into dark matter, becomes possible. As a consequence, much smaller values of A f are required in order to not overproduce dark matter.
In this case, we can use the narrow width approximation in Eq. (3.4) to find the contribution of this process to the relic density. In the case of a purely axial Z , as in Fig. 3, we have with r i ≡ m i /m Z . As we can see, when A 2 f A 2 χ , as in Fig. 3, the freeze-in contour in the Z resonance region is mostly independent of A f . For the parameters chosen in Fig. 3, the correct relic density in the Z resonance regime is found for m Z ∼ 2 GeV (m Z ∼ 200 GeV) and m Z ∼ 7.3 × 10 3 GeV (m Z ∼ 7.3 × 10 5 GeV) when m χ = 1 GeV (m χ = 100 GeV).
Note that the resonance is able to dominate the freeze-in contours even for m Z m χ , in contrast to what happens in the freeze-out contours, with narrow resonances centered at m Z ∼ 2m χ . This is because in the freeze-out case we integrate the annihilation cross-sections (= γf f →χχ /(n eq χ ) 2 , see Eq. (3.4)) up to the freeze-out temperature (T f ≈ m χ /30). In the case of freeze-in, however, we integrate the production cross-sections (= γf f →χχ /n 2 f ) up to the reheating temperature (T R m χ ).
f annihilation in the heavy Z regime When m Z m χ , m f , the s-channel happens as a four-fermion interaction, in the heavy Z regime. In the limit m 2 Z s m 2 f , m 2 χ , with m Z < T R , the contribution of the heavy regime is found by integrating the production rate over temperature up to m Z . We find If m Z > T R , the contribution of the heavy regime to the relic density would instead depend on m χ T 3 R /m 4 Z . In Fig. 3, the heavy Z regime dominates the freeze-in contour for m Z > 10 4 GeV (m Z > 10 6 GeV) for m χ = 1 GeV (m χ = 100 GeV).

Constraints on parameter space
In this section, we discuss the most stringent constraints on our parameter space. As we will see in the next section, they provide complementary bounds on the Z mass and couplings.
Indirect detection searches for dark matter could also pose limits on our parameter space. However, since they are only sensitive to annihilation cross-sections near the thermal region, they would only constrain our WIMP scenario. Such constraints would not be competitive with the experimental constraints on Z 's discussed below and are not considered here.
New species thermalized with the SM plasma at temperatures in the MeV scale can change the predictions of Big Bang Nucleosynthesis (BBN) [45,94]. In this work, we restrict ourselves to the case of dark matter candidates at and above the GeV scale. Therefore, one should only be concerned about the lower bound on Z masses. Interactions like e + e − ↔ νν, through the exchange of a light enough Z , can delay the neutrino decoupling. While a detailed analysis of such effects is beyond the scope of this work and would not change our main conclusions, we adopt the conservative bound of m Z > 10 MeV.

Direct detection
Direct detection experiments aim to identify the nuclear or electronic responses produced by the collisions between DM and the detector's target nuclei, being able to place stringent constraints and/or rule out DM models.
In our model, the scattering off nuclei takes place through t-channel exchanges of a Z . When m Z √ 2m N E R , with m N the nucleus mass and E R the recoil energy, the usual approximation of a short-range interaction via heavy mediators does not hold. In order to consider direct detection bounds with a sub-GeV Z , we use the recasted limits provided by the micrOMEGAs package [46]. We will show in Fig. 4 the constraints from XENON1T [95] on the spin-independent DM-nucleon scattering cross section, which provides the most sensitive current direct detection limits on our parameter space.

Experimental Constraints on Z Parameters
There are constraints on Z properties from a broad range of existing experiments ranging from low energy atomic parity violation measurements [48] to high energy searches at the LHC [49]. In what follows, we briefly discuss the most stringent ones, referring the interested reader to the existing literature for details. Also, our list of measurements is not exhaustive as we do not include constraints that are less restrictive than the ones we describe below. In addition, these bounds only apply to Z 's that couple to leptons and therefore do not apply to leptophobic Z 's. Our results are summarized in Fig. 4.
LHC Z bosons can be produced via Drell-Yan production, pp → l + l − X, where X represents the beam fragment jets, in hadron colliders [96][97][98][99][100][101], so that constraints can be put on Z parameters by comparing the predicted Z production cross sections for specific final states, σ(pp → Z ) × BR(Z → l + l − ), to experimental limits on these cross sections [100,101]. Both CMS [50] and ATLAS [49], and prior to this CDF [102,103] and D0 [104], have obtained such limits for specific Z models. In this work, we use the 95% confidence level experimental limits on the cross section to dilepton final states given by the ATLAS collaboration [49] based on LHC Run 2 at √ s = 13 TeV with total integrated luminosity of L = 139fb −1 . To calculate the theoretical predictions for the cross sections we use the expressions given in Ref. [105], the LHAPDF set C10 parton distribution functions [106,107], and include the 1-loop Kfactors to account for NLO QCD corrections [108,109]. NLO QCD and electroweak radiative corrections were included in the width calculations [110]. To obtain a limit on the Z couplings for a given Z mass, m Z , we take the ATLAS limit on σ(pp → Z ) × BR(Z → l + l − ) and vary the couplings for the given m Z until we obtain agreement between the predicted value and the ATLAS limit. We checked the reliability of our calculations by comparing our results with the limits on m Z for some of the models in Ref. [49]. The resulting excluded parameter space is shown in Fig. 4.
e + e − with LEP II data One can put constraints on Z 's by looking for deviations from SM predictions due to the interference with the Z in e + e − → ff . The LEP experiments, ALEPH, DELPHI, L3 and OPAL, have summarized their measurements for 130 GeV ≤ √ s ≤ 207 GeV in Ref. [51]. They give results for σ(e + e − → µ + µ − ), σ(e + e − → τ + τ − ), σ(e + e − → hadrons), A F B (µ + µ−), and A F B (τ + τ − ), where A F B are forward-backward asymmetries. We use the expressions given in [100] to calculate the predicted values for these observables. To reduce the theoretical uncertainties we use the ratio of the observable with the Z divided by the SM prediction and compare the deviation from 1 to the experimental error. We construct a χ 2 summing over all the measurements given in Ref. [51] for the observables and energy range given above to find the 95% C.L. limit on m Z . The resulting excluded parameter space is shown in Fig. 4.
BaBar from e + e − → γZ The BaBar experiment has placed limits on dark photon properties from the process e + e − → γA , followed by A → e + e − , µ + µ − , where A is the dark photon [52]. Expressions for this process are given in Ref. [111]. However, instead of attempting to properly take into account experimental details, such as detector acceptances and efficiencies, we follow a more pragmatic approach by digitizing the BaBar results [52] and rescaling them to obtain the excluded region in Fig. 4. To obtain these limits we recalculated the expression for e + e − → γZ , which resembled the expression for e + e − → γA with the substitution 2 e 2 → (V 2 f + A 2 f ), where is the kinetic mixing term between the dark photon and the SM photon and e is the electric charge. We note that BaBar reported a more recent result in e + e − → γA , where the A decays to DM, so with invisible decay products [53]. Nevertheless, the limits from this latter process are more stringent in only a few small regions of the parameter space and, therefore, we only show limits for the case with µ + µ − in the final state [52].
LHCb from A → µ + µ − The LHCb experiment has placed limits on dark photon properties from the search for dark photons produced in pp collisions at √ s = 13 TeV which subsequently decay to µ + µ − pairs [112,113]. Analogous to the previous paragraph, we rescale the LHCb limits using the substitution 2 e 2 → (V 2 f + A 2 f ). LHCb considers two scenarios, a prompt-like A search and a long-lived A → µ + µ − . Obtaining limits from the prompt-like A search is straightforward as lifetime dependent systematic effects cancel [112,113] and a Z with universal couplings to all SM particles will not alter this situation. For the long-lived A results we follow Ref. [114] and use the supplementary data from Ref. [112] to take into account lifetime dependent detector efficiencies. The long-lived A case does not give rise to any additional constraints.
Atomic Parity Violation in 133 55 Cs Atomic parity violation (APV) measurements offer some of the most precise tests of the SM electroweak theory and has been used to constrain various types of new physics including Z 's [30,48,64,[115][116][117][118][119] (see Ref. [120] for a recent review). To constrain Z 's using precision measurements of APV in 133 55 Cs we use the expressions given in Ref. [48] and the measurements given in the PDG [121]. It is straightforward to find the maximum allowed value of the Z couplings for a given value of m Z at 95% C.L.. It is important to keep in mind that the weak charge Q W is proportional to where Z is the number of protons, N is the number of neutrons and e, u and d refer to the electron and up and down quarks respectively so that we cannot obtain any constraint for the case of purely vector or purely axial couplings to SM fermions. The resulting limits are shown in Fig. 4.
(g − 2) µ and (g − 2) e The anomalous magnetic moments of the electron and muon can constrain new physics via loop contributions [64,[124][125][126], in particular, due to Z 's [90,127]. We use the expressions given in Ref. [90] to calculate the contribution to (g − 2) µ(e) from our universal Z . For (g − 2) µ , we compare this to the deviation between the average of the recent Fermi National Accelerator Laboratory muon (g − 2) measurement [128] and the Brookhaven National Laboratory experiment E821 measurement [54,55], and the SM prediction [129] (see also the electroweak review in the PDG for details [121]). The experimental average is larger than the SM prediction so we constrain the Z parameters to give agreement with the experimental average at 95% C.L.. The fitted parameter values (shown in Fig. 4 as the green bands) are completely ruled out by the neutrino-electron scattering and APV bounds. Specific charge assignments for our Z portal [127,[130][131][132] or different new physics is therefore needed to explain the (g − 2) µ anomaly. For (g − 2) e , as noted in Refs. [64,84], there is a subtlety in that (g − 2) e is the most precise measurement used to determine the fine structure constant and, consequently, the best bound to determine (g − 2) e comes from the next most precise experiment that measures α, and not from the errors from the (g − 2) e experiments themselves. Following [64,84], we use δ(g − 2) e < 1.59 × 10 −10 to constrain the Z parameters. The resulting limits are shown in Fig. 4.

Electron beam dump
We can constrain Z parameters using limits from electron beam dump experiments [64,133]. In these experiments, Z 's are produced via a bremsstrahlunglike process (e − N → e − N Z ) where the electron beam is stopped in the target with the beam products stopped by shielding. A detector looks for the decay products downstream. Thus, by comparing the expected event rate to the experimental limits, we can constrain the Z properties. To obtain limits, the couplings need to have "Goldilocks" values, i.e., not too small and not too big. On the one hand, it needs to be large enough for the Z 's to be produced in sufficient quantity but, on the other hand, if it is too large, the Z will decay too quickly for the decay products to escape the shielding. In addition, if the couplings are too small, the Z will decay beyond the detector. We follow Ref. [64], which uses the thick target approximation of Ref. [133]. One difference to note between a dark photon (A ) and our Z is that our Z can decay to neutrinos, implying that the decay width will be larger than that of an A with similar mass and couplings. The limits from the electron beam dump experiments SLAC E137 [134] and SLAC E141 [135] are shown in Fig. 4. We also considered limits from the Fermilab experiment E774 [136] but they are weaker than the BBN bound and are not shown on these plots.

Unitarity bounds
Our simplified dark matter model can violate perturbative unitarity at high energies.
In the high-energy limit, the self-annihilationχχ →χχ through a (longitudinal) Z exchange happens independently of the vector coupling [91]. The partial-wave unitarity condition in this case implies a lower bound on the Z mass: An analogous relation holds forf f →f f , with top quark self-annihilation providing the strongest limit on the (m Z , A f ) plane. However, such a limit is not more stringent than the experimental ones and in Fig. 4 we will only show the consequence of Eq. (4.1).
Unitarity bounds are of course of most concern in the freeze-out regime of our model, in which A χ is sizable. As a consequence, allowing for axial couplings can render WIMP models in tension with DM overproduction. Moreover, partial wave unitarity can establish weaker but almost model-independent upper limits on WIMP masses [137,138].
Self-annihilations of χ into a longitudinal Z violate unitarity for √ s > πm 2 Z /m χ /A 2 χ [91], such that new particles must be introduced for the consistency of the model -potentially impacting both the relic density calculation and the phenomenology. It is interesting to notice that FIMP models are usually safe from the unitarity perspective and, most importantly, the requirement of smaller axial couplings would not overproduce FIMPs. Also, for the tiny couplings involved in the freeze-in regime, the Z and χ can safely be taken to be much lighter than the particles providing their masses while restoring unitarity, which is not usually the case for simplified WIMP models.

Results
In this section, we present and discuss our results. In Fig. 4, we show the constraints on our model, as described in Section 4, for different combinations of A χ/f and V χ/f and DM masses, in the plane of Z mass and SM-Z couplings. We show the contours providing the observed dark matter relic density produced via the freeze-out and freeze-in mechanisms (solid black curves) for a universal Z portal with purely axial couplings (top left panel), purely vector couplings (top right panel) and both axial and vector couplings (bottom panels). The dashed curves set the regime boundary, with DM being produced through freeze-out (freeze-in) in the 10 0 Figure 4: The relevant constraints (coloured shaded regions), each described in Section 4, on the contours consistent with the observed DM abundance [1] (solid black curves), for both freeze-out and freeze-in mechanisms. The dashed grey curve corresponds to the regime boundary, with WIMPs (FIMPs) in the region to the left (right). In the region labeled as "non-thermal Z " (in light blue), the Z was never thermalized with the SM particles. The panels on the top (bottom) right lack this regime boundary, hence in the entire plane DM is produced through freeze-in (freeze-out). The top left panel shows the scenario with purely axial couplings, for m χ = 1 GeV, whereas the top right panel illustrates the case with purely vector couplings, for m χ = 50 GeV. In the bottom panels, we present the scenario with both axial and vector couplings, for m χ = 100 GeV, where A χ = V χ = 10 −10 on the left and A χ = V χ = 10 −3 on the right. region to the left (right) of the boundary. In the region labeled as "non-thermal Z ", the Z was never coupled to the SM fermions, as discussed in Section 3.
The top left panel of Fig. 4 illustrates the scenario where the Z couplings to both DM and SM fermions are purely axial, with m χ = 1 GeV and A χ = 10 −10 . This would be the case for a Majorana DM candidate. As shown in Fig. 3, heavier DM would require a heavier Z in order to agree with the relic abundance constraints. As one should expect, the freeze-out regime is completely excluded for such small A χ . Direct detection bounds on the spin-dependent DM scattering off nuclei, which would apply in this case, are much weaker than the other bounds and therefore not relevant. However, independent experimental limits -from neutrino-electron scattering, the beam dump experiment E137, and from the BaBar and LHCb collaborations -are currently able to probe such an elusive frozen-in DM candidate. This remarkable result extends to DM masses in the range of 10 MeV-10 TeV (with ATLAS and LEPII being able to probe the heavy DM limit and E141, the light one) and a wide range of A χ (smaller values of A χ require contours with larger values of A f ). Note that for larger values of A χ , χ is known to be a successful WIMP candidate which is able to evade the direct detection bounds (see for instance Ref. [5]).
In the top right panel of Fig. 4, only vector couplings are assumed, with V χ = 10 −10 . As opposed to the previous case, in the absence of axial couplings the production rates are too weak to allow for thermal DM in this entire parameter space, thus only the freeze-in mechanism is able to generate the observed DM abundance in this case. We also observe that the Z is more easily decoupled from the SM fermions for the same reason 5 . Moreover, the Z Z →χχ process is no longer as large near reheating, becoming negligible for the achievement of the relic abundance in this entire parameter space. In turn, thef f →χχ process is relatively independent of m Z if m Z max[m f , m χ ] (or in other words, if the process becomes Boltzmann suppressed before the resonance) and A f = A χ = 0 (cf. Eq. (3.6)). As it is already known [14,19], spin-independent direct detection bounds on the DM-nuclei scattering provide strong constraints in this case, with freeze-in already being probed for DM masses from tens of GeV to a few TeV, provided that m Z m χ . Similarly to the case of purely axial couplings, neutrino-electron scattering and bounds from the BaBar and LHCb collaborations put stringent constraints on the freeze-in contours. For larger values of V χ , smaller values of V f are needed and the freeze-in contours would also be constrained by direct detection and beam dump experiments (similarly to the case of a Z from U (1) B−L [14]). Provided that m Z < 2m χ , the conclusions above are independent of the DM mass.
In the bottom left panel of Fig. 4, we consider both axial and vector couplings of the Z to the SM fermions and DM, with a small value for the DM couplings, V χ = A χ = 10 −10 . Both freeze-out and freeze-in mechanisms can generate the observed DM abundance in the parameter space considered. For such small couplings to the Z , χ is ruled out as a WIMP by DM direct detection (limits from XENON1T) and BBN bounds, whilst LHCb and BaBar can constrain the FIMP scenario. Once again, setting the DM mass to smaller (larger) values would translate into a shift of the DM relic abundance curves to smaller (larger) Z masses. Since both axial and vector couplings are present, APV bounds are now possible and in fact are competitive with beam dump bounds for lighter DM candidates.
Finally, in the bottom right panel of Fig. 4, we show the case where the Z couples both vectorially and axially to the SM fermions and DM, with relatively large DM couplings, V χ = A χ = 10 −3 . In such a scenario, DM was able to thermalize with the SM bath in the whole parameter space, therefore being produced via freeze-out. We can see that XENON1T limits rule out most of the parameter space providing the correct amount of WIMPs, except for the well known resonance region, where m Z ≈ 2m χ (see for instance [45]) and the less explored light Z case, whereχχ → Z Z annihilations dominate freeze-out production. As discussed in Section 4.3, for such large χ − Z couplings, the unitarity lower bound on m Z shown in Eq. (4.1) becomes able to constrain part of our parameter space, as opposed to the previous cases. As a consequence, WIMPs cannot be too heavy in this case. Also note that, even thoughχχ → Z Z annihilations make WIMPs viable in a wide range of our parameter space, they violate unitarity at (not too high) energies. This renders such a simplified model less appealing, pointing towards the need of more realistic realizations (see for instance [91]).
In summary, as evident by Fig. 4, direct detection searches, experimental constraints on Z , and unitarity and cosmological bounds can currently probe and/or exclude a significant part of our parameter space in a complementary way. We have shown that, even if Z 's have very tiny couplings to dark matter and considerably small couplings to standard fermions, they are able to provide a successful freeze-in and are not invisible to current experimental searches, even in the case where they are purely axial. By considering a purely vector Z (top right panel of Fig. 4), larger values of V χ are required in order for DM to be produced through freeze-out, as compared to the pure axial case. Since in this case the correct DM abundance is produced when A f and V f are large enough, DM direct detection experiments can currently probe the freeze-in regime of the model. Future direct detection experiments such as DARWIN [139] and XENONnT [140] are also not sensitive enough to probe freeze-in if axial couplings to DM exist. We estimate that XENONnT limits probe couplings an order of magnitude smaller than the XENON1T limits for 100 GeV dark matter, which still does not reach our freeze-in contours. Nevertheless, even though direct detection fails to probe freezein in this case, the experimental constraints on Z parameters that we have considered are able to significantly constrain freeze-in DM, while leaving viable a large part of our parameter space, particularly for smaller couplings and larger m Z . Moreover, one should note that, in our scenario, the regions below the freeze-out and above the freeze-in contours are excluded by the Planck constraint, as they would overclose the universe.

Conclusions
In this work, we studied the freeze-out and freeze-in production of a Dirac fermion dark matter candidate, χ, which interacts with the Standard Model fermions, f , via a universal Z portal. We have shown how the free parameters of our model, the vector and axial-vector couplings of Z (V f , A f , V χ , and A χ ) and the masses of χ and Z (m χ and m Z ), impact the thermalization of χ and Z in the early universe. Our main results are presented in Fig. 4.
We have discussed the role of each process (depicted in Fig. 1) for the achievement of the observed relic abundance of χ (along the solid curves in Fig. 3 and Fig. 4). We found that the t-channel can only dominate the correct abundance in the presence of axial couplings, setting a lower viable value of m Z according to Planck, while s-channels dominate above this lower limit. In the absence of axial couplings (as shown in the top right panel of Fig. 4), the relic contours are dominated by s-channels and are almost independent of m Z when m Z < 2m χ .
We explored the phenomenology of this model, considering a wide range of Z masses (from MeV up to PeV) and couplings. We considered DM direct detection bounds (XENON1T), experimental constraints on Z parameters (from colliders, neutrino-electron scattering, atomic parity violation, electron and muon anomalous magnetic moments, and beam dump experiments), as well as cosmological (BBN) and unitarity bounds, as summarized in Fig. 4. Our main result is that most of these constraints can already test freeze-in in a complementary way, with viable regions such as where m Z m χ , mostly unconstrained. As expected, if the Z has purely vector couplings (top right panel of Fig. 4), existing XENON1T limits can exclude the freeze-in DM production for m Z m χ , i.e., for a sub-GeV Z . Additionally, constraints from BaBar, LHCb, and neutrino-electron scattering measurements provide much stronger bounds in the case where V χ V f , making it possible to test heavier Z 's and lighter χ's. Larger values of V χ make electron beam dump experiments sensitive to freeze-in along with direct detection bounds. Since only s-channels are responsible for freeze-in production in this case, these conclusions hold for m χ from tens of MeV up to tens of TeV provided that m Z < 2m χ .
Weakening direct detection bounds by considering V χ = 0 (as for a Majorana dark matter candidate) is known to be one of the viable options for WIMPs in simplified models. In this work we focused instead on the FIMP regime of a purely axial Z (top left panel of Fig. 4). While direct detection bounds are indeed too weak to be relevant in this case, the other experimental bounds on A f are still very stringent and able to rule out part of the viable FIMP parameter space, especially near the resonance (m Z ∼ 2m χ ). We have therefore found that for a wide range of A χ and m χ in the range 100 MeV-100 GeV, FIMPs which interact via a purely axial Z are also currently constrained by data.
In the presence of axial couplings, the production rates are stronger compared to the vector only case, making it easier for χ and Z to have thermalized with f in the early universe. The scenario of both vector and axial couplings has more stringent constraints on both FIMPs and WIMPs, now coming from atomic parity violation. If the DM couplings are very small (bottom left panel of Fig. 4), FIMPs can be tested mainly near the resonance region (for m Z < 2m χ ), similarly to the case of pure axial Z 's. Direct detection bounds are not strong enough to probe FIMPs though, since the freeze-in contours are no longer independent of m Z for m Z m χ if axial couplings exist. Larger DM couplings (bottom right panel of Fig. 4) make WIMPs viable DM candidates in very narrow regions near the Z resonance and at lower m Z values set by the t-channel contribution. The unitarity bounds prevent WIMPs from being too heavy. Also, restoring unitarity violation due to t-channels would usually require the introduction of new states which cannot be too heavy, as opposed to the FIMP case.
In summary, the proposed model offers viable dark matter candidates whose experimental signatures can already be constrained by data from a variety of complementary search strategies, showing that part of the parameter space of both FIMPs and WIMPs mediated by a Z boson can be probed at present. Interestingly, we note that, although elusive, FIMP DM can currently be probed by a variety of experiments. This motivates further work on different realizations of our Z portal, as well as the development of even more sensitive searches for new feebly interacting particles.