FIMP dark matter mediated by massive gauge boson around the phase transition period and the gravitational waves production

We study the feebly interacting massive particle dark matter whose production processes are significantly affected by the phase evolution and the complicated thermal corrections to the vector boson. We calculate the freeze-in processes to obtain the correct dark matter relic density by enumerating all the possible $1 \leftrightarrow 2$ and $2 \leftrightarrow 2$ processes. The predicted gravitational waves emitted by the first-order phase transitions and the cosmic strings are evaluated.


II. Model description 3
III. First order phase transition and the production of the gravitational waves5 A. Finite temperature effective potential 6 B. Bubble nucleation temperature T n and the percolation temperature T p 9 C. Gravitational waves from the first-order phase transition References 40

I. INTRODUCTION
Dark matter is considered to account for about 84% of the matter in our universe [1]. A natural way to explain the relic abundance is to attribute the creation of the dark matter particles to the thermal plasma in the early epoch of our universe. If the dark matter particles are initially assigned to be in thermal equilibrium with the standard model (SM) sectors, and then "freeze-out" from the plasma as the universe expands and cools down, such kind of dark matter is usually called the weak interaction massive particle (WIMP, for a review, see [2]). On the other hand, if the dark matter interacts feebly with the early plasma and is created (or "freeze-in") gradually from a void initial condition, such kind of dark matter can be called a feebly-interacting massive particle (FIMP) [3][4][5][6][7][8].
In the WIMP case, the dark matter freezes out roughly around the temperature T f ∼ mχ 26 , where m χ is the mass of the dark matter particle. Practical evaluation shows that usually the thermal corrections become negligible so that the usual methodology based upon the zero-temperature theories becomes a very good approximation to calculate the relic abundance. On the other hand, FIMPs are created much earlier when T mχ ∼ 0.3-3 in a relatively higher temperature, so the evolution of the phases at the early universe, including the phase transition processes may affect the dark matter production processes. In particular, the external leg's effective masses may vary significantly as both the vacuum expectation values (VEVs) and the thermal corrections evolve. The internal mediator's propagators for the scattering process can also change during this period. As a result, the production processes of the dark matter particles can be switched on and/or off due to the threshold effects, finally altering the interaction rates evidently during the freeze-in processes.
If vector boson participates in the freeze-in interactions [9][10][11][12][13][14][15][16][17][18], more complexities arise due to the separation of the longitudinal and the transverse degrees of freedom. For the massive vector boson, part of the goldstone boson also decouple from the longitudinal polarization degree of freedom. One has to evaluate the thermally averaged interaction rates carefully by considering all these degrees of freedom separately [19]. Furthermore, the thermally corrected dispersion relations of the vector bosons are no longer Lorentz-invariant, thus accumulates the complexity of the phase space integration.
In this paper, we rely on a model including the fermionic dark matter χ and a U (1) gauge boson A . Two Higgs singlets Φ s and Φ w are introduced to break the gauge symmetry spontaneously. As the temperature drops, Φ w and Φ s acquire VEVs successively, and phase transitions then occur to endow a mass term to the A through the VEVs. The dark matter is feebly charged under the U (1) group, so that the gauge boson participates the freeze-in processes. Φ w is also assigned with a proper charge to formulate a Yukawa coupling with the dark matter to account for the interaction between the longitudinal polarization of the A and the dark matter particle. In this paper, we consider both the thermal masses and the VEV-induced masses of the A , and their effects on the dark matter production rate.
The discovery of the gravitational waves by LIGO [20] initiates a new method to verify some new physics models through observing their predicted production of stochastic grav-itational waves in the early universe. First-order phase transitions associated with dark matter production  can emit gravitational waves to be probed by the future spacebased gravitational wave interferometers e.g., as LISA [50], TianQin [51][52][53], Taiji [54,55], DECIGO [56,57], and BBO [58,59]. In our paper, the spontaneously symmetry breaking of the local U (1) symmetry can yield the formation of cosmic strings for the scenario of Φ w Φ s , which is essential for a considerable or dominate contribution of the dark matter production rate from the longitudinal vector boson [109]. These cosmic strings then collide and self-interact to form loops, and the loops finally disappears, with the legacy of significant gravitational waves formulated via cusp, kink and kink-kink collisions [60][61][62]. We estimate the possibility to probe these gravitational waves emitted from the cosmic strings and first-order phase transitions.
This paper is organized as follows: our model is described in Sec. II; the methodology for the evaluation of phase transition and gravitational waves are given in Sec. III; we go into details for our calculation of the freeze in production process of the fermionic dark matter in Sec. IV; we comment on the phenomenological constraints on the model in Sec. V; numerical results of the dark matter and gravitational wave productions for some benchmarks are given in Sec.VI; the Sec. VII is devoted to the summary and future prospect; and some details for the readers are given in Appendix A,B,C,D.

II. MODEL DESCRIPTION
In this paper, besides the SM particles, we introduce two dark Higgs singlets of Φ s , Φ w , and one Dirac fermionic χ field that are all charged under the U (1) dark gauge group with the corresponding gauge field A µ . χ contains two Weyl components, which should always appear in pair to elude the anomaly. We impose a Z 2 symmetry under which χ is Z 2 -odd, while all the other particles are Z 2 -even. The U (1) dark charge carried by χ is denoted by t χ , and the U (1) dark charge carried by Φ s , Φ w is denoted by t s = 1 and t w . For the purpose of the freeze-in scenario, t χ 1 so χ could not directly interact with the Φ s . Assigning t w = 2t χ gives rise to the possible tenuous Yukawa coupling between χ and Φ w . After Φ w acquires the VEV, the two Weyl components of χ split, and the lighter one becomes the dark matter candidate, with its stability guaranteed by the Z 2 symmetry. With the above setups, the total Lagrangian corresponding to the dark sector is written below, where Here, , and g D is the dark gauge coupling constant. For simplicity, we define g χ = t χ g D , and g w = t w g D = 2t χ g D , so The potential term is given by where H is the SM Higgs doublet. We expand the scalar fields around their classical backgrounds as, where h, G 0 , G + , φ s , φ sη , φ w , φ wη are background fields, and the correspondingh,G 0 ,G + , φ s ,φ sη ,φ w ,φ wη are "particles". The Z 2 -odd χ can be decomposed into two Weyl spinors and the mass term can be written by where δm = √ 2y χ φ w . Diagonalizing (6) with gives rise to It would be more convenient to define two 4-dimensional Majorana spinors and then we have Therefore, the Yukawa and gauge interactions can be reduced to We are interested in the phase transition and dark matter freeze-in production process mainly around TeV-scale, and we discuss the phase transition evaluations in two scenarios. In the scenario I, v w ≈ v s ∼ TeV scale and both Φ s and Φ w appear in the phase evaluation processes. In this scenario, δm m χ so that χ 1,2 can be treated as a pair of pseudo-Dirac particles. The scenario II to be discussed is v w v s , as well as µ 2 w µ 2 s , the cosmic string produced after the φ w acquires VEV and the spontaneously breakdown of the U (1) dark , in which Φ w decouple from our TeV-scale phase transition evaluations, except its Goldstone remains φ wη which contributes to the longitudinal polarization of A . In this scenario, φ w changes little in the TeV-scale temperature, then it can be regarded as a constant, and is assigned with its zero-temperature value v w . It is then convenient to write Φ w into the nonlinear form Φ w = v w e iφwη/vw . The Yukawa term then becomes We see clearly the χ-splitting mass term above, as well as the higher order 3-and 4-point effective vertices. The v w in this scenario can become extraordinary large, thus amplify the δm to split χ 1,2 into completely two majorana fermions. Sometimes δm > m χ to induce an additional minus sign in the first eigenvalue of (8). We are going to illustrate our manipulation of it in our later discussions. The last thing we want to emphasize is that there is also a tiny coupling between χ and φ s induced by the faint mixing (denoted by V sw ) between the φ w and φ s sectors. We just parametrize such an interaction with the effective coupling

GRAVITATIONAL WAVES
In this section, we write down the methodology for calculations of phase transition and gravitational waves produced during first-order phase transition process and from the cosmic strings decay.

A. Finite temperature effective potential
For the study of the phase transition in scenario I, with the standard methodology, we adopt the thermal one-loop effective potential [63], The V 0 (h, φ s , φ w ) and V CW (h, φ s , φ w ) are the tree-level potential and the 1-loop Coleman-Weinberg potential, with V c.t 1 (h, φ s , φ w ) to keep the zero temperature vacuum from shifting. The finite temperature correction is described by the term of V T 1 (h, φ s , φ w , T ), and the Daisy-correction term V daisy Rotating the fields to expand along the φ sη = φ wη = 0 hyper plane, we obtain the tree-level potential, At the zero temperature, considering the stationary point conditions, we get In this paper, we assign v h , v s and v w as well as all the other coupling constants as our input parameters, and utilize Eq. (17) to evaluate µ 2 0,s,w . The Coleman-Weinberg contribution is given by [64] where F = 0 (1) for bosons (fermions), Λ is the MS renormalization scale, g i = {1, 1, 1, 1, 1, 1, 2, 6, 3, −12} for the {h, η, φ s , φ sη , φ w , φ wη , G ± , W, Z, T } in this model, and C i = 5/6 for gauge bosons and C i = 3/2 for scalar fields and fermions. Λ is a renormalization scale to be fixed to Λ = 3 TeV in this paper. The field-dependent Higgs mass matrix is given by, The field dependent dark photon mass is given by The slight shift of the tree-level VEVs induced by V CW is canceled by the "counter-terms" (CT) [65] with the relevant coefficients determined by, The logarithmic IR divergences encountered in (22) take the form [66][67][68][69] where φ i can be any scalar field, and G is one Goldstone mass term. We follow Ref. [65] to replace the Nambu-Goldstone boson masses with Λ IR in (22). In this paper, we adopt Λ IR = 200 GeV. The one-loop finite temperature corrections are given by [70] where the functions J B,F are with y ≡ m 2 i (h, φ s )/T 2 , and the upper (lower) sign corresponds to bosonic (fermionic) contributions, respectively. The thermal integrals J B,F given by Eq. (25) can be expressed as an infinite sum of modified Bessel functions of the second kind K n (x) with n = 2 [71], The daisy term V daisy 1 (h, φ s , φ w , T ) is given by [72,73] V daisy where the finite temperature corrections are calculated as c w (T ) = 1 12 where g 1 and g 2 are the SM U (1) Y × SU (2) L gauge couplings. The definitions of the m 2 i (h, φ s , φ w ) + c i (T ) in the mixing situation are the eigenvalues of (33), with the diagonal elements added with the c i (T ) defined in (33). The details of the mixture of the vector bosons are illustrated in Appendix. A.
For the studies of phase transition in the scenario II, the evaluation of temperature dependent effective potential is actually similar with Eqs.(14)- (33) in scenario I, with all of the φ w , λ w,sw,wh terms removed. More explicitly, after integrating out the φ w , the φ w mediated processes also converges into point-like interactions. This eliminates the φ w involved terms, while shifting the λ h,s,hs and µ 2 0,s in (15) intoλ h,s,hs andμ 2 0,s . For simplicity, we neglect the "tilde" without confusion to write down the potential from the aspect of effective theory, Therefore, the third row and column in Eq. (35) also disappears, For the field dependent mass m A , (20) becomes Notice that although g w g D , the extremely large v w φ s might still contribute significantly to the gauge boson's mass.
Since the minimum of the effective potential (h, φ s , φ w ) evolve as the temperature drops, it is necessary to study the phase evolution and transition structures of the system. We utilize both the CosmoTransitions [74] and PhaseTracer [75] by making independent programs to find out the phases as well as the transition processes among them for the cross-validation, and will only adopt the data when the results from both programs are consistent.
B. Bubble nucleation temperature T n and the percolation temperature T p For a study on first-order phase transition, one has to compute the bubble nucleation temperature T n , and the percolation temperature T p , that are usually somewhat lower than the critical temperature T c when two vacua are degenerate. The bubble nucleation temperature T n can be estimated by [76] which means that at temperatures lower than the critical temperature, at least one bubble should be created inside per-unit Hubble volume at the bubble nucleation temperature T n . The bubble nucleation rate Γ is defined by [77] where A is an O(1) constant, and S = min{S 4 , S 3 /T } is the action of the bubble solution. Usually in our model around T c,n,p , S 3 /T < S 4 so that we only display the S 3 definition where φ i (r) = (h, φ s , φ i ) is the "bounce solution" acquired from the equations of motion with the boundary conditions between the two phases φ iphase 1 and φ iphase 2 during the transition. The definition of the percolation time t p is given by [78][79][80] where Here a(t ) is the scale factor of the Friedmann-Robertson-Walker (FRW) metric, r(t, t ) is the comoving radius of a bubble given by where v b is the velocity of the bubble wall. After t p is evaluated, one can solve the T p through the equation by replacing H with 1/(2t p ) during the radiation dominant epoch. Here g is the effective degrees of freedom of the plasma, which is approximated by g 106.

C. Gravitational waves from the first-order phase transition
To evaluate the gravitational wave spectrum emitted during the first-order phase transition, one has to acquire the phase transition strength parameter of α, and the phase transition duration parameter of β, which are defined to be where ρ rad = π 2 g * T 4 /30 is the plasma energy density, and [81] is the released vacuum energy during the phase transition. Both α and β can be calculated at the phase transition temperature of either T * = T n or T * = T p for slightly different results. In this paper, we adopt T * = T p , however we still use the symbol T * , as well as the H * = H(T * ) in our following displayed equations for the purpose of generality. We then follow Ref. [77] to evaluate the gravitational wave from the first-order phase transition by summing up the contributions from bubble collision, sound wave and turbulence,

The bubble walls collision contributions
The bubble walls collision term Ω co from the "envelope approximation" results is given by [82][83][84] [110] For Jouguet detonations, we adopt the Chapman-Jouguet condition of the wall velocity v b as below [85], with the peak frequency f co locating at

The sound wave contributions
The sound wave contribution Ω sw is given by with the peak frequency being [86][87][88]: In Eq. (52), the τ sw shows the duration of the sound wave from the phase transition [89], which is calculated as where The κ ν factor in (55) indicates the latent heat transferred into the kinetic energy of plasma, which is given by [92] κ ν =

The turbulence contributions
The magnetic hydrodynamic turbulence term Ω turb is given by with the peak frequency [93] The efficiency factor ε ≈ 0.1, redshift of the frequency is obtained as

D. Cosmic strings and gravitational waves
In the scenario II, cosmic strings start to form after φ w begins to acquire its VEV. These strings collide and self-interact into loops, and then shrinks to leave us the gravitational waves formulated via cusp, kink and kink-kink collisions. The spectrum can be expressed as Particularlly, Ref. [62] transforms Eq. (60) from the Nambu-Goto string Ref. [94] into Ref. [62,95] where n = 1, 2, ..., labels the radiation frequencies ω n = 2πn/( /2). The dimensionless parameter Gµ is where v w is the VEV of the U (1) dark scalar field which spontaneously breaks down. P n is the corresponding average loop power spectrum with its numerical results adopted from Ref. [95], and C n is given by where the integration of the time parameter t has been transformed to the redshift parameter z. To evaluate the integration in Eq. (63), we need the cosmic time t(z) and Hubble constant H(z) to be expressed from the redshift parameter z to become with the current abundances of the radiation, matter and dark energy given by [1] Ω r = 9.1476 × 10 −5 , Ω m = 0.308, The function G(z) is given by 0.83, 10 9 < z < 2 × 10 12 ; while the cosmic string number density for the loops produced in radiation dominated era but survive until matter domination is given by [96]: where Γ = 50 [95], and t eq = 2.25 × 10 36 GeV −1 which is the matter-radiation equality time.

IV. FREEZE-IN PROCESSES
The dark matter particles can be produced through both decay and annihilation processes. To calculate the relic density, we rely on the Boltzmann equation. Assuming that all the other particles except the χ 1,2 are in equilibrium with the plasma, and ignoring the feedback of the dark matter particles annihilating into the plasma due to the extreme smallness of the χ 1,2 abundances compared with their equilibrium values, the Boltzmann equation is given by is the total dark matter particle number density normalize by the entropy density, γ tot is the summation over all the "rates", and x = |mχ 1 | T is the dimensionless parameter measuring the evolution of time.
To calculate the γ tot , we need to sum over all of the 1 ↔ 2 and 2 ↔ 2 processes taking into account the thermal corrections on the external legs. The thermal corrections to the dark sector particles χ 1,2 are neglected. Besides the VEV-dependent mass terms, all the other particles receive thermal corrections on their dispersion relations. These cause the complicated threshold effects, and the production rates change significantly during the freezein processes.
After the electroweak phase transition when H acquires a nonzero VEV, there will be intricate mixings of both the gauge bosons and the Higgs bosons between the SM and dark sectors. This is extremely hard to manipulate. Fortunately, in our interested parameter space when m χ 100 GeV, the freeze-in processes basically cease when x = mχ T 3, which is still well above the electroweak phase transition temperature T ew ∼ 100 GeV. Therefore, we neglect the dark matter production below the electroweak phase transition.
For the gauge boson and the SM fermions, we adopt the hard thermal loop (HTL) results to evaluate the phase space of the final states. Goldstone equivalence gauge is also utilized for the convenience to decompose the degrees of freedom of the vector boson A . The details of the HTL corrections to the vector boson and the SM fermions are illustrated in Appendix. B.
The Higgs boson masses are extracted from Eq. (14). Since we only consider the processes above the electroweak phase transition so that h = 0, and there the mixings between the SM Higgs doublet and the φ s,w vanish. Therefore, Since before the electroweak phase transition, all the elements of the SM Higgs doublet are degenerate, so m H is the mass for all of the SM Higgs bosons. Diagonalizing (70) gives two mass eigenstates mixed from theφ s,w . We use φ 1,2 to represent the two eigenstates, and m 1,2 to denote the corresponding masses. The mixing matrix elements are assigned with so that When φ s = 0 and φ w = 0, the ∂ 2 V ef f ∂φs∂φw = 0, so we need to diagonalize Eq.( 70) and calculate the mass eigenvalues and mixing matrix. In this case, the masses of φ sη and φ wη also vanish. This is because besides the gauged U (1) dark group, there is an additional global U (1) symmetry which is also broken to generate another Goldstone boson. The two massless states recombine into or one can warp the coefficients by parametrizing them into U (A ,G)(s,w) , Here φ A η connects with the A , and will be partly eaten by the longitudinal polarization of A , as will be illustrated in Appendix B. φ Gη is the Goldstone boson corresponding to the global U (1) symmetry. In the zero temperature, it might cause some phenomenology problems such as the Higgs invisible decay. We will discuss these problems later in this paper.
The assignment of the φ A η , φ Gη masses depends on the VEV structures of φ s,w . When φ s,w are both nonzero, both φ A η,Gη are massless, and if one or both of the φ s,w become zero, interactions Symbol Coupling constants The counterpart of the non-zero VEV of φ s,w is massless, while the counterpart of the zero VEV shares the same mass with the corresponding φ s or φ w . For further usage, we list all the couplings in Tab. I. Now we are ready to calculate the interaction rates γ tot . This involves evaluating the diagrams in Fig.1, 2-7 and 8 for all the 1 ↔ 2, non-SM product 2 ↔ 2 and SM product 2 ↔ 2 processes respectively. For the brevity of this section, we leave the detailed formulas in the Appendix D.

V. PHENOMENOLOGICAL CONSTRAINTS
The FIMP dark matter interacts so faintly with the SM sector, so it typically lies far beyond the ability of all the direct and indirect detection experiments. For our model  discussed in this paper, the dark Higgs sector and the A are constrained by the experimental data.
In this model, there is always a massless Goldstone boson φ Gη . The relic of such a particle behaves as a dark radiation, and can contribute to the effective number of neutrino species N eff . The SM neutrino-to-photon density ratio is given by   If the dark radiation decouples with the thermal bath at some temperature T dr when the effective degree of freedom is g * (T dec ), its density respective to the photon density is calculated to be ρ dr ρ γ = N dr 2 g * (T dec ) 4 3 .
(76) Figure 6: φ A η,Gη φ 1,2 ↔ χ 1 χ 2 diagrams. The Goldstone degree of freedom of A φ 1,2 ↔ χ 1 χ 2 can also be indicated.  If φ Gη is the dark radiation, N dr = 1, and its correction to N eff is Since φ Gη only interacts with the Higgs sector and the dark photon sector, which are massive and disappear from the thermal plasma far above the temperature T > 10 GeV, compelling φ Gη to decouple with the plasma when g * (T dec ) ≈ 100. Therefore, This is fairly safe within the Planck data [1,97,98]. The exotic Higgs bosons introduced in this model mix with the SM Higgs boson, and might alter the SM-like Higgs boson's phenomenology. According to Ref. [99], when the mixing angle between the SM-like Higgs boson and the exotic Higgs boson X sin θ hX 0.2, the SM-like Higgs boson can fit all the experimental data safely, so we set this criterion during our scanning processes. Another stringent bound in the Higgs sector is the h SM → φ Gη φ Gη , which cannot be prohibited kinematically because of the massless φ Gη . The width of the SM Higgs boson decaying into an exotic massless scalar boson X can be estimated as λ hX is the effective h-h-X-X coupling. We just estimate the bounds [97,[100][101][102]. This requires λ hX 0.02, which should be compatible with λ s,wh in order of magnitude, so we will set λ s,wh < 0.01 during our scanning process to avoid this constraint.
Finally, the A in our model might be produced at the LHC through its mixing term (D48) with the SM Z/γ bosons. Ref. [97] had summarized the detector bounds on such kind of vector bosons. For example, Ref.
[103] constrained the [σ·B]Z [σ·B]Z < 10 −5 ∼ 10 −8 depending on different Z masses. However, compared with the Z production channel, (D48) at least introduces a 2 suppression factor, not to mention the suppression due to the large A masses we adopted in this paper. As we will illustrate, we are going to fix = 10 −4 in each of our benchmark points listed in Tab. II, which is free from the Z bounds.

VI. NUMERICAL RESULTS
For all the 2 ↔ 2 processes described in Sec. IV, there exists at least one s-channel diagram. Sometimes the s-channel mediator becomes on-shell when it is above the threshold to open up a corresponding 1 ↔ 2 process. The rigorous manipulation requires to resum the one-loop "string series" of the mediator's self-energy diagrams. This modifies its thermal propagators to avoid the infinity, just similar to the familiar manipulation with the Breig-Wignar propagators in the zero-temperature case. However, due to the invalidation of the Lorentz invariance in the finite temperature, the "imaginary parts" of the s-mediators are no longer a constant for us to evaluate conveniently. Since in the freeze-in case, if any 1 ↔ 2 processes appears to be non-zero, the "off-shell" part of the 2 ↔ 2 processes are expected sub-dominant due to the extra couplings, and the"on-shell" part of the corresponding smediator should also be attributed to the 1 ↔ 2 processes, thus should be removed to refrain the double-counting. Therefore, an 2 ↔ 2 process is counted only when its corresponding 1 ↔ 2 processes disappear.
Here the λ sh , λ wh < 0.01 criterion are enough to confine the SM-like Higgs boson within the phenomenological constraints, and bounded from below criterion can be checked numerically by both the CosmoTransitions and PhaseTracer. We plot our results in Fig. 9. Among them we adopt two benchmark points for each scenario, which are called BP S1 1, BP S1 2, BP S2 1, BP S2 2, with their location in the α-β plain marked in Fig. 9, and their parameter values assigned according to Tab. II.
In scenario II, v w v s , and the spontaneously symmetry breaking around the v w -scale can create the cosmic strings, ensuing the gravitational waves relics. Notice that around the v s -scale, the Φ w sectors had been integrated out, and the only remained contribution is the y w v w within the A mass term from (36). Since y w and v w do not separate in all of the processes during the v s -scale phase transitions, one can regard them as a whole. Therefore, in the first Values Parameters BP S1 1 1 BP S1 1 2 BP S1 1 3 BP S1 2 1 BP S1 2 2 BP S1 2 3  the v s -scale phase transition induced gravitational waves still remain unchanged as g χ varies, because the product y w v w has been fixed for each of them.
For further comparison, we also plot the TeV-scale phase evolutions for all of our benchmark points in Fig. 11  h(T)>*H9@ Figure 11: Phase evolution figures for BP S1 1 (top left), BP S1 2 (top right), BP S1 2 (bottom left), and BP S1 2 (bottom right). Each of the phase evolution figure is marked with the percolation temperature T p . The BP S1 1 and BP S1 2 (top right) undergo two-step phase transitions, so T p1 , However, the cosmic strings induced gravitational wave spectrum shift as v w moves. In this paper, we adopt three sub-benchmark points BP SX Y (1,2,3) for each of the BP SX Y,corresponding to the y χ = 0.1g χ , y χ = 1g χ , and y χ = 10g χ conditions respectively. The exact values assigned to them are selected so that the dark matter relic density Ω DM ≈ 0.12 [1], and are displayed in the second table of Tab. II. (We will explain later why BP S2 1 2 and BP S2 2 3 are missing, and why the y χ is replaced with 3g χ rather than 10g χ in BP S2 1 2.) In the scenario I, this does not affect the predicted gravitational wave spectrum due to the ignorable cosmic string contributions, however in Scenario II, this will alter the values of v w to change the Gµ in (62). Therefore, we show the predicted gravitational wave spectrum for each of these sub-benchmark points as well as the sensitivities of the proposed gravitational wave experiments in Fig. 10.
With the percolation temperature T p obtained before, we can calculate the evolution of the number density of the dark matter particles. We adopt the approximation that the phase changes all in the sudden before and after t p . In scenario I, δm = y χ v w is tiny compared with m χ , so m χ is a very good approximation of the dark matter mass. However, in scenario II, evident δm = y χ v w could even cause m χ 1 = m χ − δm < 0. This contradicts with the common acknowledgement that masses are always positive. However, for the fermions, the standard manipulation is to rotate the phase of χ 1 in Eq. (8) to eliminate this minus sign in the mass terms, with the price of the accumulated evaluation complexity by casting these additional phases to the coupling constants listed in Tab. I. in this paper, we adopt a simpler but equivalent method to keep all the minus signs of the fermionic mass parameters intact in all of our evaluations. It is easy to prove that this does not affect the final results, since in the phase space integrations, m χ 1 always appears in its squared values m 2 χ 1 to eliminate the minus sign, and in the squared amplitudes, the minus sign in m χ 1 is equivalent with collecting up all the additional phase factors that makes m χ 1 positive.
In this paper, due to our limited computational resources, we start our calculation when x = 0.05. This is sufficient because the relic density evolution is not so sensitive to the starting point of x if it is sufficiently small, since the critical freeze-in temperature usually ranges within 0.1 x 3. Typically when x 5, the freeze-in processes gradually cease, and this is also confirmed by our practical calculations. Therefore, we require to iterate (68) until at least when x > 5. However, as we have stated, we find it difficult to manipulate the intricate mixings between the various Goldstone bosons after the electroweak phase transition, so we terminate our calculation before T = 200 GeV to elude the electroweak phase transition. The consistence of these two conditions removes the BP S2 2 3, in which we estimated that its corresponding m χ 1 ≈ −400 GeV, and also the BP S2 1 2, with its corresponding m χ 1 ≈ −540 GeV.
In Fig. 12 and Fig. 14, we plot the evolutions of the dark matter relic abundance for each of the sub benchmark points in Scenario I and Scenario II respectively. Fig. 13 and Fig. 15 are the corresponding evolution of the dark matter generation rates contributing to the right-handed side of (68). We list the meanings of the channel abbreviations of the γ in Tab. III. In Fig. 13 and Fig. 15 one can clearly see the alternation of the decay and annihilation  allRates Summation over all contributions. channels that dominate the total γ as the temperature evolves. This is due to the complicated threshold effects as the masses of the vector and scalar boson evolves. During the first-order phase transition, the production rate might jump discontinuously. We marked each of the x p(1,2) corresponding to the percolation temperatures T p (1,2) in all panels of Fig. 13, 15 for comparing with the phase evolution displayed in Fig. 11. Practically, when we look into the details of the annihilation channels, we found that usually χ 1 χ 2 → A * → φ 1 φ A η or χ 1 χ 2 → A * → φ 1 A L dominate the annihilation channels, because of the larger A -φ 1 -φ A coupling constants.

VII. SUMMARY AND FUTURE PROSPECT
In this paper, we construct a model where the fermionic dark matter particles feebly interact with the visible sector through the exotic vector boson. In both the dark matter production scenarios under study, it is difficult for the near-future experiments such as LISA, TianQin and Taiji to detect the gravitational wave signal produced by the TeV-scale phase transition. Meanwhile, in the scenario with the spontaneously breakdown of the dark U(1) symmetry by an extremely high scale v w , a significant gravitational wave signal is produced through the cosmic strings decay. This is well within the ability of these experimental proposals.
We enumerated and calculated all the possible 1 ↔ 2 and 2 ↔ 2 channels to produce the dark matter. The vector boson dispersion relations receive significant and complicated thermal corrections. In place of the rough estimation to treat them as something with a fixed mass for all momentum values, we applied sleeker methodology to separate the transverse, longitudinal and the partly revived Goldstone degree of freedom with different on-shell dispersion relations. The threshold effects induced by the evolution of thermal corrections as temperature drops also significantly alter the dark matter production rates. We calculate these influences and show the results for some benchmark points.
In our paper, we consider the relatively heavy TeV-scale dark matter, and avoid the troublesome mixture between the vector and Goldstone boson sectors for easier evaluation. This is valid in all our benchmark point selections, however for lighter FIMP dark matter with the mass 0.5 TeV, or the FIMP dark matter interacting directly with the SM SU(2) L ×U(1) Y gauge bosons, such effect might be inevitable. We will concern this situation and do further study on it in the future.

Appendix A: Thermal masses for SM gauge bosons
When evaluating the effective potential, one has to turn to the imaginary time formalism. The transverse and longitudinal degrees of freedom receives different corrections [72], where the script L (T ) denotes the longitudinal (transversal) mode, and W , B denote the SU (2) L triplet and hyper-charge gauge bosons respectively. Summed with the VEV-induced terms, them together, the longitudinally polarized W boson's correction is To determine the masses of the longitudinal Z and A one should diagonalize the matrix and adopt the eigenvalues to substitute the m 2 i (h, φ s , φ w ) + c i (T ) terms in Eq. (27).

Appendix B: Propagators and on-shell behaviors of the vector boson and the SM fermions
In this paper, we adopt the Goldstone Equivalence Gauge to calculate the interaction rates. In this gauge, the polarization vector is extended from 4-dimension to 5-dimension by adding up a Goldstone degree of freedom. Here we use M, N . . . to denote the 5-dimensional indices. The previous four numbers M, N, · · · = 0, 1, 2, 3 indicate the usual Lorentz time and space indices, and the last M, N, · · · = 4 represents the Goldstone degree of freedom. The full propagator can also be extended into a 5 × 5 matrix. For a particular momentum k = (k 0 , k), let us define the corresponding 5 × 5 project matrix P T , P L , where where M ± and M L are the transverse and longitudinal polarization vectors respectively. Here and in the special frameset of k = (k 0 , 0, 0, k 3 ) along the z-axis, M ± (k) = 1 √ 2 (0, 1, ±i, 0, 0). In the zero temperature, the Goldstone equivalence gauge propagator of the A is given by In the thermal plasma, adopting the "σ = β 2 -gauge" and consider the one-loop self-energy diagrams, the full propagator is given by , where a, b, · · · = 1, 2, and Π T,L,U are the factors extracted from the one-loop self-energy corrections. In the hard thermal loop (HTL) approximation, they are calculated to be where Expanded around the minima of the effective potential, Π U (0) = 0 should be satisfied. Therefore, we can use Π U (0) to estimate Π U (k), so that Π U (k) = 0. From (B8) one can immediately write down the shifted on-shell equation of both transverse and longitudinal A , Each external leg of both transverse and longitudinal vector bosons should be multiplied with the "renormalization factor" Z T,L (k), with their definitions given by A fraction of Goldstone degree of freedom has also been spewed out by the longitudinal polarization vector boson. Calculating D full,44 0 (k) from (B8), one acquires Besides the k 2 − m 2 A − Π L (k) pole corresponding to the longitudinal polarization vector boson, (B13) includes a branch cut along k 0 ∈ (− k, k) axis. The imaginary part of ∆ F GS (k) peaks near k 0 = ± k, so as an approximation we can regard these two peaks as two massless Goldstone fractions. Their "renormalization factor" is given by where and γ = c A (T ) For the leptons and quarks, we only consider the processes above the electroweak phase transition critical temperature, so the thermal corrected dispersion relation of a fermion is given by where f indicates any left-or right-handed leptons and quarks, so There are four solutions of the equations in (B16), indicating one particle, one hole, and one anti-particle as well as one anti-hole. The "renormalization factors" become Appendix C: From self-energy diagrams to 1 ↔ 2, 2 ↔ 2 interaction rates In the literature, the evaluation of the changing rate of some particle number density in the finite temperature environment relies on the imaginary part of the self-energy diagrams. In this paper, for the sake of intuitiveness we depend on the tree-level diagrams, just as everybody learned from the quantum field theory textbooks in the zero-temperature case.This appendix section aims at illustrating the equivalence between these two methods.
It is now convenient to rely on the σ = 0 gauge to calculate the imaginary part of the self-energy diagrams. Following the usual literature, let us denote 1 and 2 as the two types of vertices. The production rate of one of the dark matter particle, e.g., χ 1 is extracted from Π < (k). This is calculated from the self-energy diagrams in which the leftmost vertex is in type 2 and the rightmost vertex is in type 1. Since both the χ 1 and χ 2 are far from thermal equilibrium with the plasma, it is convenient to appoint T = 0 in all the χ 1,2 propagators only, while keeping all the other temperature terms normal in other particle propagators. Therefore,  The 1 ↔ 2 processes are extracted from the one-loop diagrams. One of the example is listed in Fig. 16. The D F 21 propagators connecting the type-1 and type-2 vertices are finally reduced to the on-shell phase-space integrals. Therefore, the one loop self-energy diagram in Fig. 16 is cut into the tree-level 1 ↔ 2 diagrams which is compatible with our previous evaluations.
The 2 ↔ 2 processes arise from the two-loop diagrams. We show two examples in Fig. 17. Cutting the propagators between the type-1 and type-2 vertices induce the 2 ↔ 2 diagrams. The left panel entails the s-channel process, and the right panel gives the t-channel processes.
Other arrangements of the vertex types are possible. In other cases, the two-loop diagrams can be cut into more than one processes. For example, Fig. 18 results in four parts if we cut all the connections between the type-1 and type-2 vertices, indicating successive real processes. However, when all the propagators connecting the type-1 and type-2 vertices become on-shell, the one-loop induced 1 ↔ 2 contributions will also arise with dominate result due to the lower perturbative orders. Therefore, we can safely omit all high-order two-loop processes when there exist non-zero one-loop contributions inducing the 1 ↔ 2 processes. This is actually what we did in our previous operations. The 1 ↔ 2 processes include the A ↔ χ 1 χ 2 and the φ 1,2 ↔ χ i χ i (i=1,2) processes. The longitudinal polarization of A involves the accessory φ A η ↔ χ 1 χ 2 terms. Before φ s,w acquire the VEVs, φ sη,wη become massive and might decay into χ 1 χ 2 . All the possible diagrams are listed in Fig. 1.
Let us calculate the A ↔ χ 1 χ 2 processes at first. In the practical evaluation, without loss of generality we rotate to the coordination that the A momentum p µ A = (E A , 0, 0, p A ). The polarization vectors are extended to 5-dimensional vectors, with the additional element assigned to the Goldstone degree of freedom, as illustrated in Appendix B. The transverse polarization vector M ± (k) = ( µ ± , 0), and µ ± (k) = (0, 1, ±i, 0). The inward longitudinal polarization vector M L (k) is given by (B4). The 4-dimensional vector part of the matrix element is given by The Goldstone part of the matrix element is Here, p χ 1 ,χ 2 are the momentum of the dark matter particles χ 1,2 which satisfy the energymomentum conservation laws p χ 1 +p χ 2 = p A . Considering the statistic and renormalization factors, the summed squared matrix elements for a particular polarization vector M t (p A ) where t = ±, L are given by where the index t does not undergo an Einstein summation. The definition of Z t is given by Eq. (B12), where T in Eq. (B12) indicates both "+" and "−". Besides the f B appeared above, we give both the definitions of the fermionic and bosonic statistic factors f F,B , The thermally averaged interaction rate can be the expressed as Notice that the p 0 A appeared in (D5) can be extracted from solving the dispersion equations where Π T,L are given by (B9), and one should adopt Π T for t = ±, Π L for t = L.
To integrate out the d 3 p χ 1 , one can boost into the frame where A is at rest. Define where θ χ 1 is the angle between p χ 1 and p A , p χ 1 A = (p 0 χ 1 A , p χ 1 A ) is the χ 1 momentum in the A 's rest frame, and θ χ 1 A is the corresponding direction after the boost. Since we ignore the thermal corrections on χ 1,2 , so their dispersion relations are not corrected, then where φ χ 1 is the azimuth angle of p χ 1 relative to the p A vector direction, and this angle remains unchanged after the boost. Substituting Eq. (D10) into Eq. (D5), we can collect all the elements to calculate the γ A ,t .
In the finite temperature case, the remained "free" Goldstone degree of freedom, which corresponds to the tachyonic branching cut in the A propagator, can be approximated to be a massless particle-like object. The direct production of χ 1 χ 2 from an "on-shell" Goldstone fraction is kinematically forbidden, so we do not need to calculate this channel. φ 1,2 can also decay into a pair of χ 1 χ 1 or χ 2 χ 2 . The dispersion relation of a scalar particle is simpler than the vector boson, and can be straightforwardly written as where m i are the masses of the two scalar eigenstates defined in (72). The matrix element is given by where i, j = 1, 2, p χ j 1,2 are the momenta of the first and second χ j particle. Then the summed squared matrix elements are given by The additional 1 2 is the interchanging factor of the identical particles.
The phase space integral can be performed with numerical algorithms, however, just as when we calculate the dark matter freeze-out processes, applying the Maxwell distribution e − p φ i T to estimate the f B ( p 0 φ i T ) can significantly simplify the phase space integral, where z φ i = m i T , K 1 is a Bessel function, and When all of the φ s,w = 0, both φ sη,wη can become massive to open the φ sη,wη → χ 1 χ 2 channel. The evaluations of the corresponding γ φsη,wη are very similar to γ φ i ,χ j , so we do not show the details here.
For the A A ↔ χ 1 χ 2 channel, the Feynmann diagrams are listed in Fig. 2 where iM µν A A ,s,χ i rs = −g µν u r (p χ i 1 )v s (p χ i 2 ) j=1,2 iM µν A A ,t,χ i rs = −u r (p χ i 1 )γ µ iM µν A A ,u,χ i rs = −u r (p χ i 1 )γ ν Here, of course p A (1,2) and p χ i (1,2) are the momentum of the external particles respectively, and they satisfy p A 1 + p A 2 = p χ i 1 + p χ i 2 . The Goldstone-vector part of the matrix element characterized by Fig. 3 is given by iM 4µ A A ,χ i rs = iM 4µ A A ,s,χ i rs + iM 4µ A A ,t,χ i rs + iM 4µ A A ,u,χ i rs , where The squared matrix elements are so the production rates are expressed as Fig. 3 and Fig. 4 also stand for the A φ Gη ↔ χ i χ i , φ A η φ Gη ↔ χ i χ i and φ Gη φ Gη ↔ χ i χ i channels. The basic steps are exactly the same with those in (D30)-(D44) to calculate the γ A φ Gη ,it , γ φ A η φ Gη ,i and γ φ Gη φ Gη ,i except that one needs to modify the momentum symbols, the coupling constants, and abolish the Z GS when a φ Gη is replacing the φ A η .
The processes φ 1,2 A ↔ χ 1 χ 2 , φ A η,Gη φ 1,2 ↔ χ 1 χ 2 and φ 1,2 φ 1,2 ↔ χ i χ i , as indicated in Fig. 5, 6 and 7 also contribute to the production rates γ φ i A ,t , γ GηA ,t and γ φ i φ j ,k , and the evaluation processes are very similar to the previous channels we discussed, except that one needs to select the appropriate couplings from Tab. I, and to manipulate the "renormalization factors" properly. In the X ↔ χ 1 χ 2 processes, the identical particle factor 1 2 is also discarded. We are not going to enumerate all of the detailed processes in this paper, however, we would like to point out that the second diagrams in both Fig. 5 and 6 are special, since the intermediate s-channel particle is the mixed propagators (B8) among φ A η and A . As an example, we calculate the corresponding matrix elements of φ 1,2 A ↔ χ 1 χ 2 to show the detailed evaluation processes.