Feebly coupled vector boson dark matter in effective theory

A model of dark matter (DM) that communicates with the Standard Model (SM) exclusively through suppressed dimension five operator is discussed. The SM is augmented with a symmetry $U(1)_X \otimes Z_2$, where $U(1)_X$ is gauged and broken spontaneously by a very heavy decoupled scalar. The massive $U(1)_X$ vector boson ($X^\mu$) is stabilized being odd under unbroken $Z_2$ and therefore may contribute as the DM component of the universe. Dark sector field strength tensor $X^{\mu\nu}$ couples to the SM hypercharge tensor $B^{\mu\nu}$ via the presence of a heavier $Z_2$ odd real scalar $\Phi$, i.e. $1/\Lambda \; X^{\mu\nu}B_{\mu\nu}\Phi$, with $\Lambda$ being a scale of new physics. The freeze-in production of the vector boson dark matter feebly coupled to the SM is advocated in this analysis. Limitations of the so-called UV freeze-in mechanism that emerge when the maximum reheat temperature $T_\text{RH}$ drops down close to the scale of DM mass are discussed. The parameter space of the model consistent with the observed DM abundance is determined. The model easily and naturally avoids both direct and indirect DM searches. Possibility for detection at the Large Hadron Collider (LHC) is also considered. A Stueckelberg formulation of the model is derived.


Introduction
The existence of Dark Matter (DM) is motivated from different astrophysical observations like galaxy rotation curves [1][2][3], bullet cluster [4], gravitational lensing [5], and cosmological observations like anisotropies in Cosmic Microwave Background (CMB) [6](for a review, see, for example [7,8]). However, we still do not know what DM actually is. DM as a fundamental particle has to be electromagnetic charge neutral and stable at the scale of universe life time.
From satellite experiments like WMAP and PLANCK [9][10][11][12][13], that measure anisotropies in CMB, we learn that DM constitutes almost 85% of the total matter content and 26.4% of the total energy budget of the universe, often expressed in terms of relic density, which provides an important constraint to abide by. Since no Standard Model (SM) particle resembles the properties of a DM particle, many possibilities beyond the SM (BSM) have been formulated to explain the particle nature of the DM, as scalar, fermion or a vector boson stabilized by an additional symmetry G DM . Amongst different possibilities, the most popular one assumes DM to be present in thermal bath in early universe due to non-negligible coupling with the SM, which eventually freezes out to provide correct thermal relic as universe expands and cools down. Weekly Interacting Massive Particles (WIMP) belongs to such thermal relic category and is widely studied due to its phenomenological richness [14][15][16][17]. However, it is also viable to assume that DM is very weakly coupled to visible sector and therefore does not equilibrate to hot soup of SM particles in the early universe and gets produced via decay or annihilation of particles already in equilibrium. Such non-thermal DM production halts after the temperature of the bath drops smaller than DM mass and the yield freezes in to provide correct relic density, see for example, [18]. DM particles which freezes in are often called feebly interacting massive particle (FIMP) and easily evades the bounds from non-observation of DM in direct or collider searches. Such a DM is mainly studied in the analysis presented here.
Effective DM-SM operators provide a model-independent framework to probe DM characteristics like relic density, direct search and collider search prospects. , often assumed to be Z 2 ) and are singlets under G SM . A heavy mediator is assumed to couple to both dark and visible sector weakly and the operators are expected to vanish when the mass of the heavy mediator goes to infinity following decoupling theorem. A complete set of such operators have been written upto dimension six assuming G SM to be SM gauge group [63,64] as well as assuming G SM ∼ U (1) EM after spontaneous electroweak symmetry breaking [65] keeping dark symmetry intact. Detailed phenomenological analysis assuming the DM to freeze-out have been carried out including the collider search prospects at Large Hadron Collider (LHC) [65][66][67][68][69][70][71][72][73].
Here, we elaborate a model where the dark sector is coupled to visible sector only via effective dimension five operator. We choose the simplest extension of the SM by Abelian U (1) X gauge group. The U (1) X vector boson is electromagnetic charge neutral and must be stable for becoming DM. The stability is guaranteed by imposing an additional Z 2 symmetry under which the dark vector boson is odd, then the kinetic mixing X µν B µν is forbidden. However a direct connection between DM and the visible sector (SM) still could be introduced if an extra real scalar (Φ) odd under the stabilizing symmetry is present. Then an operator of mass dimension five, X µν B µν Φ/Λ, is allowed. For dimensional reasons the interaction must be suppressed by an unknown new physics (NP) scale Λ. This operator has been listed in [64] and a WIMP phenomenology has recently been performed in [74]. It is worthy to mention here, even without X µν B µν Φ term, dark sector can couple to the SM, via the mixing of scalar boson (call it S) that breaks U (1) X and the Higgs doublet (H) via a portal term |S| 2 |H| 2 [35]. Here however, we will assume that the scalar S is super heavy and decouples. In addition a quartic portal interaction of the scalar Φ, Φ 2 |H| 2 is also allowed by the symmetry. The coupling is relevant for Φ being in thermal equilibrium with the SM, however fails to produce vector boson DM without the dimension five operator. It is important to note that in absence of the dimension five term, Φ becomes a stable DM candidate together with X, while the latter is completely decoupled from the SM in the limit of heavy S. With the presence of the higher dimension interaction term, X becomes stable DM, given m Φ > m X , as we assume here. We will show in sec. 4.1, that even large portal coupling of Φ 2 |H| 2 fails to contribute significantly to DM (X) production, compared to the Φ decay after Electroweak symmetry breaking (EWSB). For the consistency of Effective Field Theory (EFT), the NP scale also requires to be larger than the maximum reheat temperature Λ > T RH . Together, it is more appealing to assume that the VDM is feebly connected to the SM and it freezes-in. The paper analyzes such possibility in details. We also demonstrate the limitation of UV freeze-in which is advocated in context of effective operators [50]. We show when the reheat temperature comes closer to the DM mass scale (m) involved in production process with T RH m, massive kinematics plays an important role and IR aspects are becoming relevant.
It is worth noticing that owing to feeble DM-SM interaction to account for correct relic density in FIMP like models, the possibility of detecting such DM candidates at direct or collider searches is limited. However, if one has an extended dark sector, like we have Φ having same Z 2 symmetry as of VDM (X µ ), there can still be a possibility. We comment on seeing mono-X (where X stands for jet, photon, W, Z or H) plus missing energy signature in this framework at the upcoming run of Large Hadron Collider (LHC).
Finally, let's comment on the so called small scale cosmological problems. Even though comparison of the standard cosmological model, i.e. the ΛCDM model, with observations is very successful on scales larger than galaxies, the model has some difficulties at sub-galaxy scales, predicting, via computer simulations, too many dwarf galaxies ("missing-satellites problem") and/or too much dark matter ("core-cusp problem") in central regions of galaxies. Among possible solutions of this "small scale crisis" are e.g. models of strongly interacting DM, for VDM see e.g. [35]. However there exist also a simpler explanation of the crises. Namely, there has been an extensive investigation recently of the possibility that a realistic treatment of baryonic physics in simulations, such as supernovae feedback, stellar winds, etc. can eliminate the tension (see [75] and references therein). Therefore in this work the issue of small scale problems has not been addressed.
The paper is arranged as follows. Sec. 1 contains an introduction to VDM models. In Sec. 2 the model considered here is described and its Stueckelberg formulation specified. Sec. 3 discusses properties of the Boltzmann equation relevant for the DM production. Sec. 4 contains our findings for the DM abundance via the freeze-in and shows regions of the parameter space consistent with the observed DM abundance. In Sec. 5 we comment on experimental constraints and collider signatures of the model. Sec. 6 shows summary and conclusions. In Appendices A-D we collect useful formulae.

The Model
The minimal VDM model contains a U (1) X gauge boson denoted here by X µ . In order to enable direct interactions between X µ and the SM one also requires presence of a real scalar Φ. Both of them should be odd under a Z 2 which stabilises DM candidate, i.e. the vector boson (m X < m Φ ). Therefore the symmetry group of the model is In order to generate a mass for the dark gauge boson we also introduce a complex scalar S charged under U (1) X , which acquires a vacuum expectation value to break U (1) X spontaneously. The Z 2 transformation acts on these fields as follows: The quantum numbers under SU (3) c ×SU (2) L ×U (1) Y ×Z 2 of the new fields are tabulated in Tab. 1.
With these fields and the charges, we can write the renormalizable SU (3) c × SU (2) L × U (1) Y × Z 2 invariant scalar potential as: The total renormalizable Lagrangian then reads: where H is the SU (2) L SM Higgs doublet and D SM µ is the SM covariant derivative. The X µ field tensor and corresponding covariant derivative are defined as where g X denotes U (1) X gauge coupling. Below we are going to investigate limit of the above model when the mass of one of the physical scalars contained in the spectrum becomes very large. We expect to reproduce a version of the Stueckelberg model coupled to the extra scalar Φ and the SM Higgs doublet H. The goal is to determine, among the degrees of freedom of the considered model, the Stueckelberg scalar introduced in order to make the Stueckelberg Lagrangian gauge symmetric. Some subtleties of the limiting procedure will be addressed.

Positivity criteria
In order to formulate conditions for asymptotic positivity (for large field strengths) of the potential in Eq. (2.3) we shall first write down the matrix of quartic couplings in the basis: Now, a scalar potential biquadratic in fields is bounded from below if the matrix W is copositive [76]. Thus, the vacuum stability conditions for the potential in Eq. (2.3) are given by the Sylvester criteria for the co-positivity of W [76,77]: We emphasize that these are necessary and sufficient conditions for vacuum stability [76].

Minimization conditions and spontaneous symmetry breaking
We parametrize the scalar fields as follows: The extrema conditions for the potential in Eq. (2.3) read Hereafter we assume µ 2 H , µ 2 S , µ 2 Φ > 0 in order to generate proper symmetry breaking. We will require that the above conditions are satisfied by non-zero vacuum expectation values (vevs) of H = v h = 0 and S = v S = 0, while for Φ we require zero-vev; Φ = v Φ = 0.
The following relations are implied by the minimization conditions (2.12): We will therefore expand H and S around the non-zero vevs as follows where we have used the same notation for the fluctuations around the vacuum as earlier for the initial fields. In the expression above σ S is the Goldstone boson that constitutes the longitudinal component of the X µ , while the SM Goldstone bosons are φ ±,0 . Note that there is no potential for σ S . We have adopted a Cartesian parametrization for the doublet H together with a polar parametrization for the complex singlet S. The purpose was to find out the degree of freedom that corresponds to the Stueckelberg scalar, it will be discussed in details shortly. In Table 2 we list all possible extrema that satisfy (2.13) for v Φ = 0 together with corresponding values of the potential. There may exist three other extrema with v Φ = 0, however for the stability of Φ we are going to choose parameters that ensure v Φ = 0. We are going to find conditions that guarantee the solution #4 to be the global minimum. First we must make sure that v Φ = 0 is the only possible vev for Φ, for that purpose we will assume that for given quartic couplings we adjust is the only solution of (2.13).
Next, it turns out that . (2.16) As it will be seen shortly we assume 4λ H λ S − λ 2 SH > 0 in order to ensure positivity of masses squared, in addition λ S,H > 0 for the positivity of the potential, therefore hereby we have shown that the extremum #4 is the deepest one, and it must be the global minimum regardless what is the nature of solutions #1, #2 and #3 (refer to tabl. 2).
The mass matrix squared corresponding to the solution #4 for physical degrees of freedom expressed in the basis {h, s, Φ} reads: where it is clearly seen that only {h, ρ} mixes (as they get non-zero vevs) while Φ (the (3,3) element of the matrix) only receives contribution proportional to the vev of the other two fields. The eigenvalues of the mass matrix read: Hereafter we will adopt the convention that h 1 is always the 125 GeV SM-like Higgs boson discovered in 2012 at the LHC. Therefore m 1 = m ± and m 2 = m ∓ for h 1 heavier (upper sign) or lighter (lower sign) than h 2 . Hereafter we are going to consider the case of very heavy h 2 , i.e. m 2 m 1 . As it is seen from (2.18), for quartic couplings not exceeding perturbative limits ∼ 4π, heavy h 2 requires large v S , i.e. v S v h . It is clear that the presence of a minimum at the extremum #4 requires: (2.20) The first condition above together with the potential positivity condition (2.9) implies λ SH < 2 √ λ S λ H . Note also that 4λ S λ H − λ 2 SH > 0 guarantees positivity of v 2 h and v 2 S , see tabl. 2. Now we can now rotate the weak basis to get the mass basis via: where R is the Euler rotation matrix of the form: The mixing angle α is determined by the entries of the mass matrix as follows 2 : The potential (2.3) has 9 real parameters: (2.25) This reduces the number of free parameters in the theory to seven: All other useful relations of the parameters in the scalar potential have been furnished further in Appendix A.

Decoupling limit
Here we would like to explore the decoupling limit of a very heavy new scalar (h 2 ). From (2.18) it is clear that the limit m 2 → ∞ requires v S → ∞. In order to do that, it is useful to define: from where we see that large v 2 S corresponds to ∆ → 0. From (A.5) we obtain: Now we can investigate the behavior of the mixing angle for ∆ ≈ 0 + , it is easy to see that: . From now on we shall use the following set of parameters: 3 Then v 2 S , λ SH , sin 2α and mass parameters could be calculated and expanded in powers of m 2 as follows: Also note: With all these relations amongst different parameters of the scalar potential, assuming large m 2 , we are now going to construct an effective residual theory in the limit of large m 2 . Note that then sin 2α → 0, such that Therefore all we need to do is to expand the Lagrangian for the SM supplemented by S, X µ and Φ around the vacuum adopting the parametrization (2.14) and drop the h 2 ↔ ρ and rename h 1 by h. It turns out that the resulting effective Lagrangian reads: where the potential is independent of σ S and given by: where m h = m 1 . 3 We consider the case m2 > m1, so λH > λSM .
The kinetic terms in the limit m 2 → ∞ should be written after expanding around the vacuum and decoupling/removing ρ as follows: where the ellipsis contain all the interaction terms.

Stueckelberg Lagrangian in decoupling limit
One can easily notice that the effective Lagrangian (2.40) coincides with the standard form of the Stueckelberg Lagrangian invariant under the following transformation: In other words we have just proven that in the limit m 2 → ∞ the theory defined by the Lagrangian (2.3) reduces to the Stueckelberg Lagrangian. In addition our model is invariant under the Z 2 : There are various comments here in order. First, note that the Stueckelberg scalar is just the Goldstone boson σ S . To see this the polar parametrization of S adopted in (2.14) was crucial. A consequence of that was also the disappearance of the potential for σ S . Now let's define a current, j µ ≡ (m X X µ −∂ µ σ S ). Then note that the following potentially relevant term, j µ ∂ µ Φ, is invariant under (2.44) and therefore could be added to the standard Stueckelberg Lagrangian. However, it turns out that this operator could be omitted. It has been shown in sec. 5 of ref. [79] that terms ∝ j µ ∂ µ Φ could be removed from the Lagrangian by field redefinition: a shift of B and rescaling of Φ. Therefore hereafter j µ ∂ µ Φ will be ignored. Alternatively one could also follow arguments of ref. [80], where the author argues that the operator ∝ j µ ∂ µ Φ does not contribute to the S-matrix elements between physical states.
Since our model is gauge invariant, the quantization requires fixing a gauge. We adopt the following gauge fixing term The advantage of the adopted gauge fixing is that it cancels mixing between ∂ µ X µ and σ S . Expanding the Lagrangian one obtains eventually (2.46) In order to gauge-fix away σ S one must adopt the unitary gauge, which corresponds to ξ → ∞.
In the Stueckelberg formalism, in the unitary gauge, in presence of Φ, one could expect presence of a term like X µ X µ Φ 2 since it is allowed by symmetries. However, it turns out that the operator X µ X µ Φ 2 may only originate from dimension-6 term |D µ S| 2 Φ 2 and therefore must be suppressed by 1/Λ 2 . This is the only gauge invariant way to generate an operator ∝ X µ X µ Φ 2 . It explains why the operator X µ X µ Φ 2 can not appear as an unsuppressed dimension-4 operator, even though naively it could be added within the Stueckelberg strategy.
The decoupling limit of the scalar potential consistent with all the theoretical constraints is illustrated in Fig. 1. In the top left panel, we show the allowed region in the (λ H , λ SH ) plane for m 2 > m 1 where the colors varying from blue to yellow show larger v s . Top right panel shows the same but coloring is with respect to m 2 . Both these top panels show the decoupling limit at the outer edge of the allowed parameter space. In bottom left panel, we show λ SH as a function of m 2 at fixed values of λ H . Similarly the bottom right figure shows sin(2α) as a function of m 2 , the convergence to zero-mixing angle is clearly shown.
Since the mixing angle vanishes in the limit m 2 → ∞, i.e. sin(2α) ∝ v h /m 2 , therefore the dimension-4 interaction between dark vector and the SM disappears. Note also that since λ SΦ could be negative, therefore for λ SΦ < 0 one can always adjust µ 2 Φ so that the Φ mass squared remains at the weak scale even if v S grows. On the other hand, for λ SΦ > 0, to keep m Φ at the weak scale λ SΦ must behave as In addition, since we want to retain the vector boson mass, m X = g X v S , of the order of weak scale, therefore it is necessary that the gauge coupling diminishes as g Note also that, since for large m 2 the value of the potential at the extremum # 4 diverges as ∼ −m 4 2 /(16λ S ) therefore in order to avoid instability while preserving perturbativity we must limit ourself to large, but finite, values of m 2 .
Summarizing, to reach the Stueckelberg limit starting from the Lagrangian (2.3) a carefully adjusted trajectory in the parameter space must be adopted. An important consequence of approaching the m 2 → ∞ limit is the decoupling of the dark sector from the SM at dimension-4 operators by sin(2α) → 0 and decoupling of X µ from S by g X → 0. It should also be recalled that to avoid instability of the potential m 2 must be finite (although can be large).

Higher dimensional operator to connect DM and SM sectors
Note that the coupling λ HΦ , which parametrizes the quartic portal interactions Φ 2 |H| 2 , remains unsuppressed in the decoupling limit. Besides CP-violating operator X µνX µν 4 this is the only renormalizable (dim-4) communication between the dark and visible sectors in the limit. Leading corrections to this communication will be provided by dim-5 operators that are invariant under transformations from under Lorentz transformations and which are made of the scalar Φ, the vector X µ and possibly combinations of SM fields. It is easy to see that there are only two non-trivial operators that satisfy required symmetry conditions 5 : This operator has already been introduced in [64], where it was mentioned that such an operator can be generated at the tree level only via antisymmetric tensor mediators. Finally, one can then write down the complete Lagrangian as: where ellipsis denote interactions of the SM Higgs boson h with other SM components that are not relevant here. We note here that we necessarily assume Λ > m 2 , otherwise higher dimensional operators (neglected in this work) would appear in the scalar potential. We also adopt the following notations hereafter :

DM yield via freeze-in
It is clear from the proceeding section that the couplings ΦΦh 1 h 1 and ΦΦh 1 remain unsuppressed in the decoupling limit that we are exercising here. These interactions are ∝ λ HΦ , which is not suppressed. This is the only renormalizable communication between the dark and visible sectors. So, it is quite likely to assume that in the early universe Φ is abundant being in equilibrium with the SM (i.e. with h). Since Φ is the next lightest Z 2 -odd dark component its decays and annihilations may produce DM (i.e. X) non-thermally. This the mechanism (freeze-in) we will investigate hereafter. First, in this section, we will derive Boltzmann equations governing X production in the early universe. We will also discuss applicability of neglecting various masses while calculating the amplitudes for decays and annihilations. Before going into the details, we would like to clarify that in the following sections, in order to satisfy the EFT limit and also have a successful freeze-in, we will adopt the following hierarchy amongst the scales and masses involved in the model: where T RH is the reheating temperature at which inflaton decay products thermalize.

DM production via decay and annihilation processes
Since we are interested in the freeze-in production of the DM, hence we look for all such number changing processes with at least one DM particle in the final state. The processes that produce VDM, can easily be cooked up from interactions introduced in the preceding section and vertices collected in Appendix B. We classify all the DM number changing processes on the basis of their occurrences before electroweak (EW) symmetry breaking (EWSB) i.e. for thermal bath temperature T > T EW 160 GeV or after EWSB, i.e. for T < T EW . Due to the presence of the ΦXB vertex the VDM X can always be produced from the decay of the scalar Φ that maintains thermal equilibrium with the SM bath via the portal interaction λ HΦ |H| 2 Φ 2 . This decay channel, shown in Fig. 2, is always present before and after EWSB as m Φ > m X , which is anyway our prime assumption for the stability of the VDM. After EWSB, the decay occurs to Z, γ. Apart from the decay, we also have four 2 → 2 annihilation channels with one DM in the final state as shown in Fig. 3 before EWSB 6 . The processes include two t-channel and two s-channel graphs including Goldstone bosons (φ 0,± ). Note that, before EWSB all the SM states are massless and the Goldstone bosons (GB) are propagating degrees of freedom as the SU (2) L scalar has the form: (3.2) Figure 2: Decay of Φ → XB before EWSB and Φ → Xγ(Z) after EWSB, which contributes to the freeze-in production of X. The vertex factor shown here is the one for Φ → XB.
The dark sector fields {Φ, X}, on the other hand, are massive due to U (1) X breaking, which occurs at much higher energy scale. Due of the presence of totally anti-symmetric rank four Levi-Civita symbol in the interaction vertex ΦXB and the momenta dependent interaction vertices for the GB's, all the processes involving Goldstone bosons in t-channel and s-channel identically become zero at the level of amplitude squared. Therefore, all those processes with GB's drop out leaving only the Φ → X, B decay channel, along with the two 2 → 2 annihilation diagrams f Φ → f X, f f → XΦ (the top right and bottom left diagram of Fig. 3) for DM production via freeze-in before EWSB. However, as we shall see in subsequent sections, the decay before EWSB is sub-dominant as compared to the annihilation processes for large reheat temperature (T RH >> m). Once the EW symmetry is broken, the GB's are no more individual physical degrees of freedom, instead they become longitudinal polarizations of the charged and neutral SM gauge bosons with v h = 246 GeV. Also, the physical gauge bosons can be obtained in the mass basis by rotating the weak basis as: where c(s) w is the (co)sine of the Weinberg angle. Thus, after EWSB, the decay corresponds to Φ → X, γ(Z), while all 2 → 2 annihilation channels giving rise to DM final states are shown in Fig. 4. Due to massive propagator contributions after EWSB, the decay becomes more relevant for the determination of the DM yield, as we will demonstrate and discuss later. The decay widths and squares of the annihilation processes appearing before and after EWSB are collected in Appendix D.

Boltzmann equations for DM production
The key for freeze-in DM production is to assume that DM was not present in the early universe. In case it is produced via a decay as Φ → BX, the Boltzmann equation (BEQ) for the number density of X can be written as: where dΠ j = d 3 p j 2E j (2π) 3 are Lorentz invariant phase space elements, and f i is the phase space density of the i th particle with corresponding number density being: where g i is the number of internal DOFs. It is important to note that we assume negligible abundance of X as the initial condition, also we disregard Pauli-blocking/stimulated emission effects, i.e. we assume 1 ± f i ≈ 1. Indeed it has been assumed that Φ's are in equilibrium with the thermal bath (SM). Similarly, the BEQ for DM production via generic annihilation process i, j → X, k (with one DM in the final state) reads [14]: The BEQ in Eq. (3.6) can be rewritten as an integral over the CM energy as [18,81]: where 2 in the limit m a,b → 0. Next we define the yield Y X ≡ n X /s, as a ratio of DM number density n X and the comoving entropy density in the visible sector s. The BEQ corresponding to the decay in terms of the yield Y X can be written in the differential form as: where we have defined: as the decay width of Φ. It is possible to express Eq. (3.8) in terms of the dimensionless quantity x ≡ m X /T and the reaction density γ D for decay as: where γ D , called reaction density, is defined in Appendix C.
For the case of annihilation one can similarly write: Note that the sum over i, j, k indicates all the possibilities of producing DM following Figs. 3, and 4. Again, in terms of the reaction density defined in Appendix C, one can express the yield due to annihilation as: Following Eq. (3.8) and Eq. (3.11) the total yield due to decay and due to annihilation can be written as: (3.13) The maximum temperature available to the process is what we call reheat temperature T RH . Also, we note that the processes after EWSB, are different from those before. Therefore, taken all such processes together, the yield at temperature T 0 can finally be written as: where the first parenthesis corresponds to contribution before EWSB while the second one describes the after EWSB production and M (a)bEWSB stands for the amplitude for processes appearing (after) before EWSB. Note that for the annihilation processes we have considered the massless approximation which makes the expression less complicated. One can, equivalently, express the BEQ in a more general manner in terms of the reaction densities utilising Eq. (3.10) and Eq. (3.12) as: This is rather more common and convenient way of parametrization. In Sec. 4 we will be using these reaction densities to compare DM yield before and after the EWSB.

UV limit and limitations
In this section, we demonstrate the difference between massless and massive limit of DM production cross-section and therefore we will be able comment on limitations of the Ultra Violet (UV) freeze-in advocated in [50,82]. Here we limit ourself to the period before EWSB so all the SM particles are massless. The masses of the dark sector Φ and X are assumed to be of the same order, typically m Φ ∼ m X ∼ m ∼ O(1TeV). Hereafter the "massless limit" refers to zero-mass approximations, i.e. both SM and dark masses are zero. Our task in this section is to estimate size of mass effects, i.e. contributions to the yield that depend on the dark masses. Since we limit ourself to the temperatures above T EW therefore Eq. (3.14) can be simplified 7 as: Here onward we shall refer to M bEWSB as M. 7 The contribution from Φ decays are negligible here.
For strictly massless case the integration over s can be performed analytically, so the result reads where θ min is an angular cutoff necessary to avoid singularities that appear in the forward direction for t-channel diagrams, in the following θ min = 10 −2 will be adopted. Now, we would like to estimate the difference between the massless and massive limit. For that, let us define:  Note that ∆Y (T EW ) is linear in m as it obviously should vanish in the limit m → 0. We should also note that Eq. (3.23) is process independent and can be written in this particular form whenever the matrix element squared can be expressed in the form of Eq. (3.18). Now, let us find out the condition under which ∆Y becomes as large as Y massless , as that will dictate the condition for which massless limit can be no longer adopted for the UV freeze-in scenario. This can simply be found out by equating Eq. (3.23) and Eq. (3.17):

(3.24)
This implies that, as long as the masses involved in the freeze-in process are approximately less that three times of the reheat temperature, one can overlook the masses and the yield can be computed in the massless limit. This, in other words, justifies the fact that for large reheat temperature the massless limit is a good approximation for obtaining the UV freeze-in yield. In Fig. 5, we demonstrate the difference between the yield obtained for massive and massless case of DM production for all the processes before EWSB put together. We plot Y X (T = T EW ) as a function of T RH and see that massive case sharply differs from the massless case (black line) at small T RH , while they exactly merge as the reheat temperature becomes large. In the bottom panel we show the same feature in terms of ∆Y Y m=0 for different choices of m X,Φ .

DM relic abundance via freeze-in
As described in details in the last section, within the freeze-in paradigm, the DM yield is controlled by the annihilation and/or decay of SM as well as dark sector particles. Following Eq. (3.15) we write down the BEQ governing the DM yield as [35]: where Y X = n X /s. H is the Hubble parameter given by H = 1.66 √ g ρ T 2 /M pl and x = m X /T is a dimensionless quantity to parametrize the temperature of the thermal bath. As mentioned earlier, the γ's denote the so-called reaction density for different particles annihilating (decaying) to the DM. The detail expressions of the reaction densities for 2 → 2 annihilations and 1 → 2 decays are given in Appendix. C. In order to compute the DM yield, we solve Eq. (4.1) with the initial condition Y X ≈ 0 at large T i.e., small x in accordance with the usual FIMP set-up 8 . By solving Eq. (4.1) one can obtain the total DM yield Y X at the present epoch i.e., Y X (T 0 ). The relic abundance of X at present temperature can then be obtained via: We must also remind that the PLANCK [13] allowed relic density reads: which we will use to find the relic density allowed parameter space of the model. Since in our case the connection between the dark and the visible sector proceeds via a nonrenormalizable interaction, the DM abundance is usually characterized by UV freeze-in [50,82] limit, where the DM abundance is sensitive to the reheat temperature T RH of the universe and NP scale Λ only. This is in sharp contrast to the Infra-Red or IR freeze-in scenario where the two sectors communicate via renormalizable operators, and the DM abundance is set by the IR physics i.e., the yield becomes maximum at low temperature, typically at T ∼ m X [18]. Now, the reheat temperature T RH is very loosely bounded. Typically, the lower bound on T RH comes from the measurement of light element abundance during Big Bang Nucleosynthesis (BBN), which requires T RH 4.7 MeV [83]. The upper bound, on the other hand, comes from (a) cosmological gravitino problem [84,85] in the context of supersymmetry, that demands T RH 10 10 GeV to prohibit gravitino over production and (b) simple inflationary scenarios that require T RH ∼ 10 16 GeV [86,87] for a successful inflation. The reheat temperature, thus, can be regarded as a free parameter for our analysis. As we have already shown in the preceding section, when reheat temperature drops close to the masses involved in the annihilation or DM production process, massive kinematics start playing a key role in the yield and UV limit can not be trusted. Therefore our analysis will be divided into two regimes, depending on the scale of the reheat temperature (T RH ): • T RH >> m i , where the reheat temperature lies way above different masses that appear in the model.
• T RH m i , where the reheat temperature lies close to the masses in the theory.
In the following sub-sections we will consider the above two cases. It is important to point out that our calculations of cross-sections have all been done analytically (see Appendices) and then the relic density is found out by solving BEQ in Mathematica numerically.

Large reheat temperature: T RH >> m
In this sub-section the reheat temperature is assumed to be much larger than the dark sector masses. We would like to mention that computing the reaction densities, we consider both 2 → 2 annihilation channels and Φ → XB decay channel before EWSB, while after EWSB we only take into account the Φ → XZ(γ) decays for reasons elaborated soon. Also note that, all the SM fields are massless before EWSB but the dark sector fields (m Φ , m X ) are massive irrespective of the era, thanks to U (1) X breaking at very high scale. First we compare reactions densities for annihilation and the decay as a function of x ≡ m X /T . Often we use a dimensionless variable r = m X m Φ < 1 to illustrate scan results. In Fig. 6 we show individual contributions of annihilations and the decay to the reaction densities as a function of x keeping the DM mass, the Φ mass, and the coupling α(α) fixed. The vertical dashed-dotted line represents EWSB x EW = m X /T EW . Here we clearly see that the reaction densities due to the s-channel and t-channel annihilation processes are much larger than that due to decay before EWSB, while after EWSB the reaction density falls to a very small value.
The suppression of the decay before EWSB can be understood comparing the analytical forms of reaction densities for annihilation (in massless limit) and decay as follows: where σ ann (s) ∝ α 2 (1 + β 2 ) × const. was assumed. It then follows that for T >> m Φ =⇒ q 1 and for T ∼ m Φ =⇒ q ∼ 1. We however, would like to caution the reader that above formula is not strictly valid at T close to EWSB when massive kinematics become important for computing annihilation cross-sections as elaborated in Sec. 3.3. This implies that the DM is dominantly produced before EWSB, while the yield accumulated after EWSB is negligibly small. This is a typical feature of UV freeze-in where the maximum yield production happens at high temperature. An evolution of the total reaction density taking into account all annihilation processes occurring before EWSB together with the decay after EWSB for different sets of {m Φ , m X } is shown in Fig. 7. It is clear from the discussion above that decay contribution before EWSB is very small, while the contribution from decay is dominant over annihilation processes after EWSB. Note here that the vertical dashed-dotted lines in different colours represent EWSB at x EW = m X /T EW corresponding to those m X values chosen for illustration. Figure 7: Total reaction density over the whole range of temperature before and after EWSB taking into account both annihilation and decay before EWSB, but only decay after EWSB. We consider all SM states to be massless before EWSB, while dark sector particles are massive. The dashed-dotted vertical lines correspond to x EW = m X T EW for each choice of DM mass. In the left panel we show densities for different choices of m Φ = {50, 100, 300, 500} GeV (in red, green, blue, black, respectively) for fixed r = m X m Φ = 0.2, while in the right one they are shown for different r = {0.2, 0.4, 0.6, 0.8} (in red, green, blue, black, respectively) for fixed m Φ = 100 GeV. In both cases we choose α = 10 −14 GeV −1 , β = 0.1.
Also note that for m Φ 100 GeV, the reaction densities in two regimes before and after EWSB can be distinctively identified by a step, while that for larger m Φ they are continuous. This is because with temperature dropping, the 2 → 2 production cross-sections drops and around T EW the decay contribution dominates over them, resulting a continuous curve before and after EWSB. However, after EWSB, the decay final state changes from Φ → X, B to Φ → X, γ(Z). Therefore when m Φ < m X + m Z , one of the decay processes, Φ → X, Z, become kinematically forbidden showing a distinct drop in the reaction density. The values of different parameters chosen for illustration are mentioned in figure inset.
The effect of large reheat temperature and large reaction densities at high temperature is reflected in the DM yield shown in the top panels of Fig. 8. In both scans we keep T RH = 10 8 GeV. In the top panel of Fig. 8 we see that the DM yield builds up quickly with x (i.e. with lowering temperature) and reaches its maximal value already at very high temperature T ∼ T RH . Then it freezes-in immediately producing an yield that remains constant till T = T 0 2.73 K. The asymptotic value of the yield, Y 0 ≡ Y (T 0 ), directly implies the PLANCK observed DM relic abundance via (4.2). As seen from (4.2) each choice of the DM mass requires appropriate Y 0 what results in the splitting of the colored curves at large x observed in Fig. 8. On the other hand each Y 0 requires the couplings α,α tuned appropriately, as shown in the legend of Fig. 8. The left top panel corresponds to m Φ = 500 GeV, while the top right panel to slightly smaller value m Φ = 100 GeV. The vertical dashed lines show the locations of EWSB, although its effect on the final yield is invisible. In the bottom panels of Fig. 8 we show contours corresponding to the central value of the PLANCK observed relic abundance (Ωh 2 0.1199) in the α −m X plane for fixed m Φ = 500 GeV and 100 GeV and two different values of T RH . The left and right lower panels correspond to the reheat temperature T RH = 10 8 and 10 6 GeV, respectively. The relic abundance is obtained following Eq. (4.2).
Since Ω X ∝ m X , we see for larger DM mass smaller α is required to compensate for the over abundance. Note, that in each panel the kinematical condition m Φ m X is obeyed. As expected, for the same DM mass, growing T RH requires lower couplings. We would finally note that to find yield in such a scenario, the masses in all reactions can be safely neglected and the processes after EWSB contributes negligibly. DM production via annihilation or decay processes can also be compared (at T 0 ∼ 0 GeV) by naive dimensional analysis as advocated in [18,88]: where T F I denotes the characteristic freeze-in temperature scale at which the yield reaches the constant value (see e.g the plateau in upper panels of Fig. 8) for DM production via annihilation process or decay, which are not quite the same. For decays, the freeze-in temperature can be assumed to be the mass of the decaying particle, i.e. T F I ∼ m Φ , which is used in the second step of the above analysis. On the other hand, for DM yield produced via annihilation process, the freeze-in temperature can be assumed to be the highest temperature available for the process, i.e. T F I ∼ T RH . Therefore, for T RH >> m Φ ensures that annihilation contribution dominates over decay contribution in the final yield for UV limit. This brings us to an asymptotic formula of the yield in UV limit, which can simply be written as: neglecting the decay contribution and is validated in Fig. 11, as we explain in a moment.

Low reheat temperature: T RH m
This case is more interesting. It turns out that when the reheat temperature drops down to ∼ TeV scale, then processes that take place after EWSB are relevant and contribute significantly to final yield. After EWSB all the SM states also become massive, and hence in order to get meaningful results, all masses shall be kept. As we shall see, in such a case the IR freeze-in starts showing up, i.e. the freeze-in takes place at a temperature T ∼ m X . Before solving BEQ, let us estimate first the hierarchy of reaction densities. For illustration, in Fig. 9, we compare reaction densities for the Φ decay and those for annihilations into X, γ(Z) final states. For the later final state two annihilation diagrams contribute: (a) t-channel annihilation h, Φ → X, Z via Z boson mediation and (b) s-channel process h, Φ → X, γ(Z) via Φ mediation. We consider all the states involved in these two processes to be massive. Note that the s-channel amplitude for h, Φ → X, γ(Z) is proportional to λ HΦ , therefore it could be amplified. We show the reaction density as a function of x for decay as the black curve in all the figures, while we choose three values of λ HΦ = {0.1, 1, 4π}, shown respectively in red, green and blue, for the annihilation processes. It is seen that for small m X 100 GeV, after EWSB the decay dominates over annihilation even when the portal coupling is close to its limiting perturbative value i.e. 4π. However for m X 400 GeV and λ HΦ 4π the s-channel annihilation starts dominating over decay for x < ∼ 7. This is expected since growing DM mass causes phase space suppression of the decay width. However, even for m X 499 GeV for large enough x again the decay dominates over annihilation. Therefore, it is fair to conclude that as long as λ HΦ O(1), one can safely ignore all the annihilation processes even for large DM mass. In the bottom panel of Fig. 9 we show variation of Φbranching ratio to X and photon or Z-boson. For light X, Φ decays into X, Z dominate while for r = m X /m Φ > ∼ 0.008 Φ decays mostly into Xγ. Now let us find relic density in the low T RH regime. In the left and right upper panels of Fig. 10 we have shown the DM yield for several values of m X as a function of x for T RH = 10 4 GeV and 10 3 GeV, respectively. In order to satisfy the relic abundance, for smaller T RH , we need larger α as Y X ∼ T RH α 2 . Note that the left panel already shows formation of IR-like behavior of the yield that is typical for low T RH : the "second slopes" that appear for larger x in the left panel (T RH = 10 4 GeV) reach their corresponding plateaus at the same location as in the right panel for T RH = 10 3 GeV. In the right panel, the typical UV-like sudden DM production disappeared leaving only after-EWSB IR-production (decays) at much larger x. The IR behavior is, of course, more prominent for smaller T RH . In the lower panels we have shown curves in the m X − α plane that reproduce proper DM abundance for a fixed T RH . Comparing with the high T RH regime, Fig. 8, one can observe that since Y X ∼ T RH α 2 the required α had to be two orders of magnitude smaller than here.
Finally, let us present an approximate formula for the yield for the case of low reheat temperature (T RH ∼ m). Proceeding in a similar way as in the previous subsection, we can estimate the contribution to Yield from annihilation and decay via dimensional argument as in Eq. 4.5. However, we need to remind that now the freeze-in temperatures (T F I ) are roughly the same for both annihilation and decay when T RH ∼ m Φ , resulting in similar decay and annihilation contributions to the yield, i.e. Therefore, the final yield for such a situation can be written as: where T F I ∼ m characterizes dark sector mass.

Summary results
Effects of varying reheat temperature for DM yield evolution has been shown in the top panels of Fig. 11 for two different sets of dark sector masses. For T RH = 10 8 GeV (shown by the red curve) i.e. for the T RH m, we observe yield that follows typical UV freeze-in pattern and becomes maximum at T ∼ T RH . With gradual decrease in T RH , although the characteristic UV freeze-in is still visible at smaller x, however the yield also builds up at larger x and final freeze-in occurs at T ∼ m X , as shown by the blue (T RH = 10 4 GeV) and black curves (T RH = 10 3 GeV). For T RH = 1 TeV the IR characteristic of freeze-in is more prominent. Note that, all the parameters chosen in these plots reproduce the observed relic abundance and hence the yields for the same m X at low x converge to the same asymptotic value, as they indeed must do according to (4.2). Colorful curves in the upper panels correspond to different values of α adjusted so that in spite of varying T RH , the asymptotic value is the same and corresponds to the observed abundance. 9 For large T RH one requires a smaller α to obtain the right abundance to compensate the effect of larger integration region. We also note that the yield Y (x → 0) values in these the two upper panels are different due to different choices of DM masses. In the bottom panel of Fig. 11 we show curves in the α − T RH space that imply proper DM abundance. It is interesting to note that for T RH 1 TeV, the relic density becomes independent of the reheat temperature as the IR freeze-in dominates over UV freeze-in. Beyond 1 TeV the effective coupling α must decreases with grow of T RH for a fixed DM mass in order to satisfy the central value of the PLANCK observed relic abundance as Y X ∝ T RH α 2 . Also for a fixed T RH , larger DM mass requires smaller α simply because Ω X ∝ m X following Eq. (4.2).
5 Signature of the model As we have already demonstrated, the dimensionful Φ−B−X vertex (α(α)) is constrained by relic density of DM to be α(α) 10 −12 GeV −1 for T RH 10 3 GeV. Therefore, assuming c ∼ O(1) 10 , one can conclude that freeze-in of the VDM requires the NP scale Λ > ∼ 10 12 GeV at least or higher for larger reheat temperature. If the effective operators (2.47) where n-loop generated, then we would roughly conclude that Λ > ∼ 10 12−2n GeV. Hereafter we assume that the Lagrangian (2.47) is indeed tree-level generated. Therefore, the phenomenology of the model is severely constrained. First of all, note that there is no XX − SM vertex, therefore no elastic scattering of the DM against nuclei is possible. So, this model easily avoids stringent constraints from non-observation of DM scattering at direct search experiments. DM can only scatter off nuclei inelastically with Φ in the final state as shown in left panel of Fig. 12.
Since m Φ > m X , hence such an inelastic scattering is forbidden even if the mass difference 9 The initial condition for all our solutions of the freeze-in BEQ assumes no DM at TRH, i.e. YX (mX /TRH) = 0. On the other hand α is adjusted so that the observed abundance is satisfied, i.e. for the same mX the curves in upper panels of Fig. 11 converge to the same value for x → 0. 10 As mentioned in ref. [64] the Lagrangian (2.47) can be generated at the tree-level by integrating out anti-symmetric tensor mediators, then indeed c ∼ O(1). [89]. On the other hand, due to the presence of Φ − X − γ vertex, the DM pair annihilation may give rise to monochromatic X-ray line (right panel of Fig. 12) but such photon flux will be hugely suppressed by 1/Λ 2 and can not account for the, say, galactic-center gamma ray excess as observed. Hence this model in its freeze-in realization of DM can not be probed from either direct nor indirect DM search experiments.  However, since the portal coupling λ HΦ is unconstrained by existing observations and as argued before this coupling can be as large as λ Hφ ∼ O(1), hence the non-standard scalar Φ can be produced at the collider via Higgs mediation (see Fig. 14). Once these φ's get produced at colliders they will eventually decay via Fig. 2 to DM and SM final states like Z, γ. This is the same channel that also gives rise to DM production via freeze-in. For particular choice of the effective coupling α that gives rise to right relic abundance, the Φ remain stable over the detector length. This is shown in Fig. 13, where we see that for α = 10 −14 − 10 −13 GeV −1 the average decay length of Φ is L Φ = cτ Φ 100 km, where c = 3 × 10 8 m/sec. In such a case the Φ's basically escapes LHC detector and gives rise to missing energy ( / E T ), which can be constructed out of the recoil of an initial state radiation (ISR) of a gluon, γ, W ± , Z, H as where the sum runs over all visible objects that include leptons and jets, and unclustered components. Therefore, the model can finally produce monojet 11 plus missing energy signal that has extensively been searched at the LHC [90-92] as a vanilla DM signal particularly for Higgs portal DM models. However, usually, when one produces DM that is connected with the SM via a Higgs portal, then the coupling is tightly constrained from direct search. Therefore, such signals are pretty small and submerged into huge SM background. In our case, as mentioned, the coupling (λ Hφ ) can be large and can produce a significant number of such mono-X signal events, that may be of interest for next run of LHC. It is worth noting that decays of Φ have to be completed before the onset of Big Bang Nucleosynthesis (BBN) [14,93], so that it does not alter the standard BBN picture. Therefore here we will require that τ Φ < ∼ τ BBN ∼ 1 sec which, for fixed m Φ and m X < m Φ puts a lower bound on α. This has been illustrated in Fig. 14 where we show, for m Φ = 100 GeV and 500 GeV, regions allowed by the BBN constraint (cτ = 10 5 km). Concluding one can see that usually α > ∼ 10 −15 GeV is allowed, while the region m X ∼ m Φ is forbidden for any value of α.

Summary and Conclusions
A vector boson DM weakly coupled to the visible SM sector via dimension-5 operator has been presented and the parameter space allowed by observed relic DM density has been found.
The advantage of the model is the absence of tree-level elastic DM scattering against nuclei and a double suppression of present time DM annihilation in, for instance, dwarf galaxies. Therefore this scenario easily and naturally satisfies existing experimental constraints. The model contains a dark sector composed, in a unitary gauge, of a massive vector X µ , and a real scalar Φ. X µ is a gauge boson of spontaneously broken extra U (1) X gauge symmetry. The vector X µ and the scalar Φ are odd under a Z 2 symmetry introduced to stabilize the DM candidate, X µ . Φ is assumed to be heavier than X µ . The SM sector is extended by an extra heavy neutral Higgs boson that decouples when its mass goes to infinity, as we assume here. The lowest dimensional operator responsible for DM-SM interaction are 1/Λ X µν B µν Φ and 1/ΛX µν B µν Φ. It has been shown that the model can be formulated in a Stueckelberg-like fashion as a limit of the SM extended by the U (1) X gauge symmetry together with a complex scalar charged under the U (1) X (needed to spontaneously break the symmetry) and a real scalar Φ.
We have investigated a possibility of DM production via a freeze-in mechanism through decays of Φ and annihilations including Φ. It turned out to be convenient to consider two distinct regimes of the reheat temperature. The first one is when the reheat temperature is significantly higher than masses involved in the production process. This situation mimics the case of UV freeze-in, when the production happens mostly before EWSB and all processes after EWSB are insignificant. However the situation alters, when reheat temperature (which can be thought of a free parameter, being very loosely constrained by BBN) drops to lower values close to the mass scale (m) typical for the dark sector. It has been shown that UV freeze-in, although advertised to describe the case of freeze-in production of DM in EFT formalism, is not fully correct, massive contributions start playing an important role and effects of IR freeze-in i.e. DM yield building even up to lower temperature (T ∼ m) starts showing up.
In order to predict properly the observed DM abundance, the scale of the dimension-5 operators must be large Λ ∼ 10 12 − 10 16 GeV depending on the DM mass m X , the reheat temperature T RH and an underlying mechanism for the generation of the relevant effective operators. The huge size of Λ implies that at the lowest level of perturbative expansion neither elastic scattering off nuclei is allowed nor present time annihilations of DM in e.g. centers of galaxies are possible. However, it turns out that LHC collider signals mediated by Higgs boson exchange are possible, gg → H * → ΦΦ. Since the scale of Λ required by the DM abundance is large Λ ∼ 10 12 − 10 16 GeV the heavier scalar Φ is effectively stable at the detector length scale and hence can produce mono-jet, photon, Z, W ± or H events accompanied by missing energy drifted away by pairs of Φ bosons.The signal cross-section could be quite substantial as the portal coupling between Φ and the SM remains unconstrained.
Finally, we must mention that a freeze-out possibility of the same model can also be thought of. In that case, the phenomenological signatures will become richer. In contrast to the case considered here the freeze-out scenario implies constraints that are more difficult to satisfy [74].

B Relevant vertices
Adopting the Lagrangian of the model in Eq. (2.48), one finds relevant vertices and propagators collected in the table 3. Here the notation have usual meaning, for example, g 1,2 are the gauge couplings corresponding to U (1) Y and SU (2) L gauge groups, respectively. c v and c a are defined as: c f v = T 3L −2 sin 2 θ w Q f and c f a = T 3L , where T 3L is the SU (2) L isospin quantum number and Q f is the charge of the SM fermion f concerned.

C Reaction densities
For a 2 → 2 annihilation channel the reaction density is defined as: where a, b(1, 2) are the incoming (outgoing) states and g a,b are corresponding degrees of freedom. Here f eq i ≈ exp −E i /T is the Maxwell-Boltzmann distribution. The Lorentz invarint 2-body phase space is denoted by: The amplitude squared (summed over final and averaged over initial states) is denoted by |M a,b→1,2 | 2 for a particular 2 → 2 scattering process. The lower limit of the integration over s is s min = max (m a + m b ) 2 , (m 1 + m 2 ) 2 .
For a 1 → 2 decay process the reaction density is given by: (C.2) D Expressions for squared amplitudes before EWSB

D.1 t-channel annihilation before EWSB
The spin averaged amplitude squared for f Φ → f X process is given by: where N c = 1(3) for the SM leptons (quarks). Also note that all the SM fermions are massless. In the limit m Φ = m X = 0 this reduces to a relatively simplified form: Corresponding annihilation cross-section is given by:

D.2 s-channel annihilation before EWSB
The spin averaged amplitude squared for f f → XΦ process is given by: 4) with N c = 1(3) for SM leptons (quarks). In the limit m X = m Φ = 0 this reduces to: Corresponding annihilation cross-section is given by: (D.6)

D.3.1 Before EWSB
The amplitude squared for the Φ → X, B decay is given by: The resulting decay width can be written as: