Strongly self-interacting vector dark matter via freeze-in

We study a vector dark matter (VDM) model in which the dark sector couples to the Standard Model sector via a Higgs portal. If the portal coupling is small enough the VDM can be produced via the freeze-in mechanism. It turns out that the electroweak phase transition have a substantial impact on the prediction of the VDM relic density. We further assume that the dark Higgs boson which gives the VDM mass is so light that it can induce strong VDM self-interactions and solve the small-scale structure problems of the Universe. As illustrated by the latest LUX data, the extreme smallness of the Higgs portal coupling required by the freeze-in mechanism implies that the dark matter direct detection bounds are easily satisfied. However, the model is well constrained by the indirect detections of VDM from BBN, CMB, AMS-02, and diffuse $\gamma$/X-rays. Consequently, only when the dark Higgs boson mass is at most of ${\cal O}({\rm keV})$ does there exist a parameter region which leads to a right amount of VDM relic abundance and an appropriate VDM self-scattering while satisfying all other constraints simultaneously.


I. INTRODUCTION
In spite of increasing astrophysical and cosmological evidence for the existence of the dark matter (DM) [1,2], the nature of DM remains a mystery. According to the dominant paradigm DM consists of collisionless, cold particles that successfully explain the large scale structures in our Universe. However, collisionless cold DM predictions obtained by N-body simulations face some difficulties known e.g. as the cusp-vs-core problem [3][4][5][6] or the too-bigto-fail problem [7][8][9] when confronted with precise observations at the dwarf scale. However, it has been shown that the presence of sizable DM self-interactions with σ DM /m DM = 0.1 ∼ 10 cm 2 /g has the potential to alleviate such a tension [10][11][12][13][14][15][16][17], even though the DM selfinteractions are constrained to be σ DM 1 cm 2 /g by measurements at the cluster scale [18][19][20][21][22][23].
Such large DM self-scatterings naturally arise if there is a light particle mediating the DM interaction and the corresponding cross-section is enhanced by non-perturbative effects [24][25][26][27][28][29][30][31][32][33]. One immediate consequence of this light mediator scenario is that the DM self-interaction cross section is velocity dependent [14][15][16][19][20][21][22], which allows for the signals at the dwarf scale to evade the constraints from the galaxy clusters. A simple way to realize this scenario is to introduce a model, where DM is generated via the dark freeze-out mechanism in which it predominantly annihilates into a pair of light mediators. Nevertheless, it has recently been shown in Refs. [34,35] that this secluded DM model [36] is severely constrained by the DM indirect detection. A way to avoid these problems is to consider a DM production mechanisms different from the conventional freeze-out. One possibility is the freeze-in mechanism [37,38] (see i.g. Ref. [39] for a recent review and the complete references therein). It is found in Refs. [37][38][39][40][41] that the final DM relic density is determined exclusively by the main DM production channels at the freeze-in temperature and it is not sensitive to many details of DM evolution at higher temperature, which guarantees the predictability of this mechanism.
The freeze-in as a production mechanism for self-interacting dark matter was analyzed in [42][43][44][45][46][47]. Notably, the case of light-mediator was discussed in ref. [44] within the model of Hidden Vector DM with dark SU (2) gauge symmetry [48], where it has been found that the scenario with keV mediator agrees with experimental constraints. It has been also noticed that if decays of the mediator into e + e − are allowed, its significant abundance and large lifetime cannot satisfy bounds from Big Bang Nucleosynthesis (BBN), so that this region of the parameters is excluded.
In this work, we study an abelian version of vector dark matter (VDM) models [48][49][50][51][52][53][54][55][56][57][58] in which the VDM particle with mass of O(GeV ∼ TeV) couples to the SM sector only through the Higgs portal. We take into account recent bounds from BBN, CMB and discuss possibility of constraining the model with FERMI-LAT, AMS-02, diffuse γ/X-Ray and direct detection LUX data. In the case of indirect constraints on DM annihilation, we include the effect of Sommerfeld enhancement. We also take into account consequences of electro-weak phase transition in calculation of DM production. The dark Higgs boson of the VDM model is assumed to be so light that it can induce large self-interactions to solve the small-scale structure problems. We focus on the scenario in which the VDM is produced by the freeze-in mechanism. The main question that we address is if there exist a region in the parameter space that can generate the right VDM relic abundance and appropriate DM self-scatterings while still satisfying all the direct and indirect detection constraints. After scanning over the parameter space we conclude that if the mediator h 2 is too light to decay into e + e − , then indeed all the constraints can be satisfied together with correct relic abundance and appropriate DM self-scatterings. The necessary mediator mass is of the order of O(keV).
Our results agree with those found in [44].
The paper is organized as follows. In Sec. II, the VDM model is presented. The VDM production via freeze-in is discussed in Sec. III, with a special attention to the effects of the electroweak (EW) phase transition. Then we discuss the VDM self-interactions in Sec. IV.
Sec. V and Sec. VI are devoted to constraints from DM direct and indirect detection experiments. The numerical results are presented in Sec. VII. Finally, we give a brief summary in Sec. VIII. Some useful formulae are collected in Appendix A.

II. THE MODEL
Following Refs. [53,54], we introduce a dark U (1) X gauge symmetry and a complex scalar S which is neutral under SM gauge group but has unit charge under this U (1) X symmetry.
We further assume an additional Z 2 symmetry, under which the gauge boson X µ and S transform as follows: which is just the charge conjugate symmetry in the dark sector. It forbids the kinetic mixing between the SM U (1) Y gauge boson B µ and X µ , X µν B µν , ensuring stability of X µ . Therefore, the relevant dark sector Lagrangian is given by where H is the usual SM Higgs SU (2) L doublet, and the covariant derivative of S is defined as D µ S ≡ (∂ µ + ig X X µ )S with g X being the corresponding dark gauge coupling constant.
Note that the quartic portal interaction, κ|S| 2 |H| 2 , is the only connection between the dark sector and the SM, so in the limit κ → 0 the two sectors decouple. Also, the mass term of S has the negative sign compared with the usual scalar field, so that it can induce the spontaneous symmetry breaking (SSB) of the gauge U (1) X . By minimizing the scalar potential of the model, we can obtain the vacuum expectation values of the usual SM Higgs Note that S can be always assumed real without compromising any generality, therefore the discrete symmetry (1) remains unbroken as needed for the stability of X µ .
After the SSB happens, the dark gauge boson obtains its mass m X = g X v S via the dark Higgs kinetic term, and both scalar fields can be written as By expanding the scalar potential up to the second order, the mass squared matrix M 2 of the two physical scalars (φ H , φ S ) T is given by With the following orthogonal transformation of scalars, we can define the mass eigenstates (h 1 , h 2 ) T with their masses (m h 1 , m h 2 ), where θ is the mixing angle with s θ ≡ sin θ and c θ ≡ cos θ. As a result, we have the following relations: In the freeze-in mechanism, the dark sector composed of X and h 2 never thermalizes with the visible SM sector, so that the portal interactions κ or s θ should be very tiny. As is evident from Eq. (6), the h 1 boson is mostly SM-Higgs-like, while h 2 is almost the dark Higgs φ S . We have found that the most convenient choice of input parameters which specify the models is (m X , m h 2 , κ, g X ), together with the already known parameters v H = 246 GeV and m h 1 = 125 GeV.

III. VECTOR DARK MATTER RELIC DENSITY VIA FREEZE-IN
Within the freeze-in mechanism, the standard assumption is that the initial abundances of the VDM and the dark Higgs h 2 after reheating are assumed to be negligibly small, which is possibly a result of the reheating process itself or another mechanism. Furthermore, the Higgs portal coupling should be very tiny so that the dark sector can neither thermalize itself nor be in equilibrium with the SM sector. When the VDM mass is smaller than the EW phase transition temperature T EW 160 GeV [59,60] its abundance is mainly controlled by various SM particle annihilations and/or decays that contribute to the collision term of the following Boltzmann equation where Y X = n X /s is the DM yield defined as a ratio of DM number density n X and the entropy density in the visible sector s. The parameter x ≡ m X /T describes the SM sector temperature T , H is the Hubble parameter, and γ i ≡ σv i (n i eq ) 2 denotes the so-called reaction density [41] for the SM particles annihilation into VDMs (for γ f we sum over all SM fermions). The last term represents the SM-Higgs-like h 1 decays to a VDM pair when this channel is kinematically allowed. Since here m h 2 m X , no corresponding decay term for h 2 appears. In this project, the model is implemented within LanHEP [61,62] and calculations of the cross sections and decay rates are performed adopting CalcHEP [63]. Definitions of reaction densities, obtained cross sections and decay rates are collected in Appendix A.
It is interesting to note that all of the reaction densities are proportional to the square of the Higgs portal coupling κ 2 with no dependance on g X 1 , which explains why we have decided to use κ instead of sin θ as a parameter. Also, due to the assumed mass hierarchy m h 2 m X , the value of m h 2 influences the resulting DM abundance very weakly. Hence, the prediction for VDM relic abundance depends mainly on two parameters m X and κ. Since the freeze-in mechanism is IR dominated [38,41], the VDM relic density is dictated by the h 1 → XX decay rate. We present the resulting evolution of the VDM yields Y X in Fig. 2, which illustrate typical features of the freeze-in mechanism.  However, when the VDM mass is much larger than the EW phase transition temperature T EW , the VDM abundance stops increasing before the EW phase transition. In this case, the SM gauge symmetry SU (2) L ×U (1) Y is not broken, so that only the tree-level diagram shown in Fig. 3 can generate VDM particles. Hence, the Boltzmann equation can be simplified to: argument in Ref. [38,40], the yield could be estimated as where the first relation follows from the dimensional argument with M Pl being the Planck mass. σ(T FI ) is the total cross section of the SM particle annihilation at the freeze-in temperature T FI , which is simplified to be σ ∼ κ 2 /T 2 FI . We have also used the relation T FI ∼ m X , which can be understood as follows. When m X > T EW , as it has been mentioned above only the channel HH † → XX contributes to VDM generation. It becomes ineffective as the temperature drops below m X , since then the SM Higgs doublets do not have enough kinetic energies. On the other hand, for the case with m X ≤ T EW , the VDM freeze-in process is dominated by the annihilations of particles which are lighter than the VDM.
Similarly, when the SM plasma temperature decreases below m X , the VDM yield ceases to grow any more due to the fact that these channels are no longer kinematically allowed.
Concluding, the freeze-in temperature is expected to be around the VDM mass, T FI ∼ m X , in the present scenario. Then it is easy to derive from Eq. (10) that the predicted VDM relic density Ω X h 2 ∝ Y X m X should only depend on κ whereas the dependence on m X are cancelled out, which is manifested as a flat line in Fig. 4. However, if the VDM is lighter than a half of the visible Higgs mass, the decay channel h 1 → XX dominates, so that the VDM yield should be where the decay rate should be Γ h 1 →XX ∼ κ 2 m h 1 , and the freeze-in temperature in this case is T FI ∼ m h 1 at which the density of the visible Higgs h 1 is greatly suppressed by its Boltzmann factor. Hence, the VDM relic density is Ω X h 2 ∝ κ 2 m X , which results in the Fig. 4. Finally, note that the small but abrupt rise of κ at the m X = 160 GeV represents the EW phase transition effect due to the sudden change of the main VDM production channels.
In order for the freeze-in mechanism to work, it is required that the dark sector neither thermalize by itself nor with the SM sector. It is easy to check that the portal coupling κ implied by the VDM relic density is so tiny that it is impossible for the dark sector to equilibrate with the visible one. However, the non-thermalization of the dark sector by itself is not guaranteed. When the number densities of the VDM and h 2 accumulated via freezein become large enough, it is probable that the dark sector process XX → h 2 h 2 would be cosmologically efficient, which would soon change the number densities of VDM and h 2 to form a dark plasma with a common (and in general different from the SM) temperature.
Therefore, one should ensure that thermalization in the dark sector cannot take place and the appropriate condition can be coded by the following inequality [41,44]: where σ(XX → h 2 h 2 )v , n X , and H represent the thermally averaged cross section for VDM annihilations into h 2 pairs, the number density of VDM, and Hubble parameter, respectively, all of which are evaluated at the freeze-in temperature T FI . Note that σ(XX → h 2 h 2 )v is proportional to the dark gauge coupling α 2 X , so that it should not be suppressed in the parameter space where the DM has large self-interactions. Thus, the condition in Eq. (12) is not easy to be satisfied in the present scenario and therefore it constraints substantially the freeze-in parameter space as shown below.

IV. VECTOR DARK MATTER SELF-INTERACTIONS VIA A LIGHT MEDIA-TOR
It is well known that the cosmological small-scale structure problems, such as the 'cusp vs. core' and the 'too-big-to-fail' problems could be ameliorated if DM self-interaction was sufficiently strong at the dwarf galaxy scale [10][11][12][13][14][15][16][17], the required value of the cross-section where σ T ≡ dΩ(1 − cos θ)dσ/dΩ is the so-called momentum transfer cross section between DM particles. However, DM self-scattering cross-section as large as σ T /m X 10 cm 2 /g is not allowed by observations at the cluster scale with the typical constraint σ T /m X < 1 cm 2 /g [18][19][20][21][22].
A possible strategy that may generate large DM self-interaction is to introduce a mediator which is much lighter than the DM particles. In the VDM model, the elastic DM scattering is mediated by an exchange of the two Higgs scalars, h 1 and h 2 . In the limit of small mixing, the h 1 -mediated contribution is negligible due to sin α and large h 1 mass suppression. In contrast, XXh 2 coupling is not suppressed by small mixing and, in addition, it is much lighter than the VDM particle, therefore h 2 can act as a light mediator which might be capable to amplify the self-interaction. When α X m X m h 2 with α X ≡ g 2 X /(4π) the finestructure constant in the dark sector, the perturbative Born approximation is applicable in which the dominant t-channel h 2 -exchange to the transfer cross section as follows [29]: where v is the relative velocity in the VDM two-body system. Nevertheless, beyond the Born range, h 2 is much lighter than α X m X so that the nonperturbative effects would become important and give rise to the following attractive Yukawa potential: Note that due to such nonperturbative corrections, the DM self-interactions have the nontrivial dependence on the VDM velocity. When the range of the potential characterized by 1/m h 2 is much larger than the VDM de Broglie wavelength 1/(m X v), i.e., m X v m h 2 , this part of parameter space is well known as the classical regime, for which analytic fitting formulas for σ T [28,29,31,64] are available in literature. In our numerical calculations, we adopt the more recent improved analytic expressions provided in Ref. [31]. On the other hand, if m X v m h 2 , the VDM self-scatterings can be enhanced by several orders of magnitudes due to the formation of the quasi-bound states. This region of parameter space is usually denoted by the resonant regime. In this work, we obtain σ T in this regime by closely following Ref. [29] to solve the non-relativistic Schrödinger equation with the potential in (15). Moreover, it has been found [28,29,31,64] that, with the presence of the non-perturbative effects, the VDM transfer cross section σ T is enhanced more significantly as the relative DM velocity becomes small. Such a velocity dependence of VDM self-scatterings is very appealing, since it helps the VDM model to solve the small-scale structure problems for the dwarf galaxy scale with a typical velocity v ∼ 10 km/s while evading the strong constraints from the galaxy clusters with v ∼ 1000 km/s. More recently, a more careful analysis of DM self-interactions from a light (pseudo-)scalar has been presented in Ref. [35], where a more appropriate definition of the momentum transfer cross section σ T is given and the possible correction from the u-channel light mediator exchange is investigated. However, it is seen in Ref. [35] that such corrections lead to very small modifications in final results so that we neglect them and follow the conventional formula from Refs. [28,29].

V. DIRECT DETECTION OF THE VECTOR DARK MATTER
It is usually claimed that the DM direct-detection experiments do not provide relevant constraints for models in which the DM particles are mainly produced by the freeze-in mechanism since the DM nuclear recoil cross sections are suppressed by tiny portal couplings.
However, in the present scenario, the spin-independent (SI) VDM-nucleon (XN) scatterings are mediated by the two neutral Higgs bosons h 1,2 , and thus it is possible that the crosssection is greatly enhanced by the small mass of the light mediator h 2 . This feature is clearly reflected by the corresponding formula for the differential cross section of the XN scatterings with respect to the momentum transfer squared q 2 , where v is the VDM velocity in the lab frame, µ XN ≡ m X m N /(m N + m X ) is the reduced mass of the XN system, and is the total cross section for the XN scattering with the effective nucleon coupling f N ≈ 0.3 [65][66][67]. Compared with the usual definition of the SI independent DM-nucleon cross section in the literature, Eq. (16) has an additional form factor G(q 2 ) defined as which encodes the effects of the light mediator h 2 . It is clear that, for the heavy mediator case with m 2 h 2 q 2 ∼ 4µ 2 XN v 2 , the factor G(q 2 ) will be reduced to 1, i.e., we will recover the conventional XN contact interaction, and the usual experimental constraints can be applied.
But when m 2 h 2 q 2 , the XN differential cross section in Eq. (16) will have extra q 2 dependence characterized by the G(q 2 ), thus modifying the corresponding nuclear recoil spectrum and, in turn, the final fitting results. Therefore, we need to re-analyze the experimental constraints in the latter case.
The strongest constraints on the direct detection of the VDM come from the LUX [68], PandaX-II [69] and XENON1T [70]. In the present work, we use the LUX 2016 dataset as an illustration of the SI direct detection limits to the VDM model since PandaX-II and XENON1T datasets would give the similar results. Due to the modification of the DM nuclear recoil spectrum caused by the light mediator h 2 , we follow the simplified analysis methods presented in Ref. [71,72]  understood that the h 2 mass is cancelled out in the final expression in Eq. (16) in this parameter region. However, even though it is remarkable that the LUX upper limit of κ reaches the order of 10 −10 for large VDM masses, it is not able to give a meaningful constraints on the freeze-in region of our model. Thus, in the following, we will not consider the direct detection constraints any more.

VI. INDIRECT DETECTION CONSTRAINTS ON VECTOR DARK MATTERS
Phenomenology of indirect detection of VDM crucially depends on properties of the mediator h 2 , such as its mass m h 2 , lifetime τ h 2 and dominant decay channels. Since we are interested in the light h 2 which could give rise to the large enhancement of VDM selfinteractions, we will limit ourself to m h 2 100 MeV. Thus, the parameter space is naturally divided into two regions: (i) m h 2 ≥ 2m e and (ii) m h 2 < 2m e , where m e is the electron mass. In the former region, the dominant h 2 decay channel is e + e − pairs, while only the diphoton mode is kinematically available in the latter case. Consequently, the light mediator lifetime τ h 2 is different in these two regions. Specifically, 10 4 s τ h 2 10 12 s in region (i) while τ h 2 10 12 s in region (ii), which is illustrated in Fig. 6 for a typical VDM mass m X = 100 GeV and a Higgs portal coupling κ = 2.09 × 10 −11 consistently with the DM relic density (see Fig. 4). Analyzing constraints from DM indirect searches, we will consider these two regions separately. Such a late decay of h 2 would produce e + e − pairs with sufficient energy that would spoil the predictions of abundances of various elements [73][74][75][76]. We adopt the most recent results from Ref. [76] where the authors also studied the BBN effects triggered by decays of dark Higgs bosons produced by the freeze-in mechanism. It is seen from Note that the result in Ref. [76] was obtained in the limit of κ → 0 and v S → ∞ while keeping θ fixed, so the 2 → 2 processes involving top quarks predominate the h 2 production via freeze-in. However, in our scenario, the Higgs portal coupling does not approach zero. The most important contribution to h 2 density arises from the SM-like Higgs decay h 1 → h 2 h 2 , which is more efficient than the top quark annihilations and top-gluon inelastic scatterings. Therefore, we expect that h 2 is more abundant in the our model, which leads to even stronger constraints. In other words, the application of dark Higgs results in Ref. [76] here leads to the conservative constraints. and S is the s-wave Sommerfeld enhancement factor given by [29,[83][84][85] S = π a sinh(2πac) with a ≡ v/(2α X ) and c ≡ 6α X m X /(π 2 m h 2 ). Since the velocity of the VDM was very small during the photon last scattering, we can use the value of S saturated in the vanishing velocity limit. Due to the large mass hierarchy between the VDM X and the mediator h 2 , the CMB upper limit in Fig. 8 of Ref. [86] for the one-step cascade with the e + e − final state can be applied for the VDM annihilation cross-section.
• AMS-02: The local annihilations of VDMs into h 2 pairs decaying to e + e − in the final state can lead to an excess of positron flux in cosmic rays [87][88][89]. Therefore, the absence of such an excess would give rise to a strong upper bound on the VDM annihilation cross section. Currently, the most precise measurements of the positron flux [90] and positron fraction [91] come from the AMS-02 Collaboration. By taking into the account the Sommerfeld enhancement factor in Eq. (19) with typical VDM velocity v X ∼ 10 −3 in our Galaxy, we can take the AMS-02 positron flux constraints from Ref. [86] for one-step cascading VDM annihilations. Note that the AMS-02 results are reliable only down to the DM mass ∼ 10 GeV, since the positron flux spectrum lower than 10 GeV would be affected significantly by the solar modulation so that the constraints in this range would be uncertain.
• Dwarf Limits from Fermi: The VDM annihilations in the dwarf spheroidal galaxies provide bright γ-ray sources in the Milky Way, and are thus expected to be probed and constrained by the Fermi Gamma-Ray Space Telescope [92]. In the present model with m h 2 > 2m e , most γ-rays are generated by the final-state radiation from the mediator decay h 2 → e + e − γ, which follows the VDM annihilation XX → h 2 h 2 . However, due to the suppression factor from radiative corrections compared with the dominant decay channel h 2 → e + e − , the constraints from Fermi shown in Ref. [86] are much weaker than the corresponding ones from CMB and AMS-02. Therefore, we do not show dwarf limits from Fermi in our following numerical results.
We now turn to the indirect search constraints for the VDM with the mediator mass m h 2 < 2m e , in which h 2 decays dominantly in the diphoton channel, and the lifetime is typically longer than 10 12 s. As mentioned before, for such a light h 2 , the BBN constraints can be evaded as shown in Ref. [76].
• Dwarf Limits from Fermi: Since h 2 → γγ is the dominant h 2 decay we expect that there should be strong constraints from measurements of γ-rays by Fermi Gamma-Ray Space Telescope [92]. However, note that the signal region for each dwarf is defined as the one within an angular radius of 0.5 • . For the 15 dwarfs used in the Fermi-LAT analysis, their distances from the Earth range from 32 kpc to 233 kpc. Thus, due to the fact that h 2 propagates at the speed of light without any scatterings in the range of a dwarf, the h 2 will spend, at most, the time of O(10 11 s) traveling inside the signal region from the center of the dwarf. In other words, it is too short in time for h 2 to decay inside the signal region. As a result, the Fermi-LAT constraints in Ref. [92] cannot be adopted directly in our case.
• CMB: For τ h 2 > 10 12 s, h 2 would have a large abundance at the time of recombination.
Also, the high-energy photons from h 2 decays would ionize and heat neutral hydrogen after recombination, and hence distort the CMB anisotropy spectrum. Consequently, recent measurements of the CMB by Planck [77] can provide strong constraints on h 2 properties [93,94]. We adopt the recent lower bound on the decaying DM lifetime τ 0 for the diphoton final state shown in Fig. 7 of Ref. [95] to obtain the following constraint for the VDM model: where Ω h 2 h 2 is the current h 2 relic density generated via the freeze-in mechanism if h 2 were present today without decays. In fact, h 2 might decay well before. The constraint is actually for h 2 abundance at the epoch of recombination, not today. We only use the present DM relic density as a reference to quantify the h 2 density fraction at the recombination period. Moreover, the expression on the right-hand side in Eq. (20) is just an approximation and the true formula should be Ω h 2 h 2 /(Ω h 2 h 2 +Ω X h 2 ). However, the h 2 density is always smaller than that of VDM due to the assumed mass hierarchy, so that Ω h 2 h 2 in the denominator can be neglected. Note that the exclusion limit in Ref. [95] extends to the DM mass of 10 keV, so we ignore the CMB constraints below this VDM mass in our numerical calculations.
• Diffuse γ/X-Ray Bounds: When the lifetime of h 2 is larger than the present age of the Universe τ U = 4.3 × 10 17 s, the h 2 particle contributes to the present DM relic density even though it is not absolutely stable. The only decay channel h 2 → γγ could be constrained by the accurate measurement of the diffuse γ/X-ray background.
Following Ref. [44,[96][97][98], we adopt the conservative lower limit on the h 2 lifetime as where Ω h 2 is the relic abundance of h 2 generated via the freeze-in mechanism.

VII. NUMERICAL RESULTS
Having discussed all of the VDM signals and constraints, we can put everything together to see if we can find a region in the parameter space where large DM self-scatterings for the scale of dwarf galaxies can be compatible with the VDM relic density and all of indirect search constraints. Note that there are four free parameters in our original VDM model, so that if the Higgs portal coupling κ is fixed as shown in Fig. 4 by the requirement that the VDM relic density constitute all of the DM in the Universe, we can plot the parameter space in the m X -α X plane with fixed values of m h 2 . The final results for some typical values of m h 2 in the Regions (i) and (ii) are presented in Fig. 7 and Fig. 8, respectively.   The situation changes a lot for the Region (ii) as shown in Fig. 8, since the indirect detection constraints are all imposed on the decay process h 2 → γγ, rather than VDM annihilations. Here we only consider the freeze-in region below the thermalization curve. In both panels, the signal regions for the dwarf galaxy scale are all in the Born and classical regions, in which the part with a small hidden gauge coupling α X and a light VDM corresponds to the Born region while the band with large values of α X and m X to the classical region. The discontinuities in both plots represent the mismatch of the analytical formulae around the boundary of these two regions. By detailed calculations, it is found that all of the signal regions for m h 2 10 −2 MeV are constrained tightly by observations of CMB and diffuse γ/X-rays, as illustrated by the left panel of Fig. 8. Only when the h 2 mass is reduced to O(keV) a small parameter window opens, in which the VDM mass is around O(GeV) and α X is in the range 10 −9 ∼ 10 −6 , as is seen clearly in the right panel of Fig. 8. In this appendix, we collect the formulae for the relevant SM particle annihilation cross sections as well as SM-like Higgs h 1 decay rate, which are involved in the calculation of dark matter relic density in the Universe via the freeze-in mechanism. Since the SM EW phase transition has a substantial impact on the VDM production, we consider the annihilation and decay channels in broken and symmetric phases, respectively.

EW Symmetry-Broken Phase:
Quark Annihilation to VDMs: where we use the SM Higgs boson width Γ h = 4.15 MeV [99,100] to regulate the SM-like Higgs mass pole singularity.
Lepton Annihilation to VDMs: W Bosons Annihilation to VDMs: Z Bosons Annihilation to VDMs: SM-Like Higgs Bosons h 1 Annihilaiton to VDMs: where we have only kept the leading-order terms in the double expansion of κ and m 2 h 2 /m 2 X . SM-Like Higgs Boson h 1 Decay to VDMs: When the VDM mass is smaller than a half of the the SM-like Higgs mass, the VDM can also be produced by the decay of h 1 , with the decay rate as follows: where we have used the definition of κ in Eq. (7) and the approximation that c θ ≈ 1.
SM-like Higgs Boson h 1 Decay into h 2 's: where we only keep the leading order in the expansion of κ and m h 2 /m h 1 because m h 2 m h 1 .

EW Symmetric Phase:
SM Higgs Doublet H Annihilation to VDMs where m H and m s are the masses of the SM-Higgs doublet and the dark Higgs φ S defined in Eq. (4).
Note that in our derivation of the Boltzmann equation in Eq. (8) after and before the EW phase transition, we have used the so-called reaction density γ i for various channels.
For the annihilation channels, the reaction density is defined as [101] γ(a b → 1 2) ≡ dp a dp b dp 1 dp 2 f eq a f eq b (2π) 4 δ 4 (p a + p b − p 1 − p 2 )|M(a b → 1 2)| 2 where a, b (1, 2) represent the incoming (outgoing) particles with g a,b as their respective degrees of freedom, and f eq i ≈ e −E i /T is the Maxwell-Boltzmann distribution. Here dp ≡ d 3 p/(2π) 3 (2E), |M| 2 are the amplitude squared summed over quantum numbers of the initial and final states without averaging, and s min = max[(m a + m b ) 2 , (m 1 + m 2 ) 2 ].

ACKNOWLEDGMENTS
We would like to thank Kai Schmidt-Hoberg, Laura Covi, and Chao-Qiang Geng for useful discussions. This work is supported by the National Science Centre (Poland) research project, decision DEC-2014/15/B/ST2/00108.