Relativistic freeze-in with scalar dark matter in a gauged B − L model and electroweak symmetry breaking

We explore relativistic freeze-in production of scalar dark matter in gauged B − L model, where we focus on the production of dark matter from the decay and annihilation of Standard Model (SM) and B − L Higgs bosons. We consider the Bose-Einstein (BE) and Fermi-Dirac (FD) statistics, along with the thermal mass correction of the SM Higgs boson in our analysis. We show that in addition to the SM Higgs boson, the annihilation and decay of the B − L scalar can also contribute substantially to the dark matter relic density. Potential effects of electroweak symmetry breaking (EWSB) and thermal mass correction in BE framework enhance the dark matter relic substantially as it freezes-in near EWSB temperature via scalar annihilation. However, such effects are not so prominent when the dark matter freezes-in at a later epoch than EWSB, dominantly by decay of scalars. The results of this analysis are rather generic, and applicable to other similar scenarios.


Introduction
The null-results from a number of dark matter direct detection experiments motivate to explore alternate dark matter production mechanisms. One of the most well-motivated dark matter production mechanisms is freeze-in [1] production of dark matter. In this framework, the dark matter is feebly coupled with the Standard Model (SM) particles, in general with particles in equilibrium and thereby referred as feebly interacting massive particle (FIMP). Due to very suppressed interaction, the FIMP dark matter never attains thermal equilibrium with particles with which they are feebly coupled. The suppressed interaction further gives natural explanation for the non-observation of any direct detection signal. The dark matter in freeze-in scenario is produced from the decay and/or annihilation of SM and beyond Standard Model (BSM) particles which are either in equilibrium [2][3][4][5][6][7][8][9] or also JHEP05(2021)150 freezing-in along side the dark matter [10]. We explore freeze-in production of dark matter in extended gauged B − L model, where we address few of the subtlety of the production.
The gauged B − L model [11][12][13] is one of the most appealing, yet minimal theory descriptions, that explain small SM neutrino masses. The model includes three right handed neutrinos (RH neutrinos) required for anomaly cancellation, one B − L gauge boson, and a complex scalar field. The scalar field acquires vacuum expectation value, and breaks the B − L gauge symmetry. The B − L gauge boson, as well as, the heavy neutrinos acquire their masses due to B − L symmetry breaking. The light neutrinos, on the other hand, acquire their masses via seesaw [14,15], with their masses inversely proportional to the B − L symmetry breaking scale. A scalar particle with B − L charge can be accommodated in this model which serves as the dark matter candidate by suitable choice of B − L charge. The freeze-out scenario for this model has been explored in [16][17][18][19][20][21], along with other phenomenological implications. The late decay of RH neutrinos are explored in scalar dark matter freeze-out scenario [22,23]. In this context the fermionic dark matter has also been explored [24][25][26][27][28][29][30]. The freeze-in scenario along with neutrino mass and leptogenesis has been studied for this model in [9]. The freeze-out scenario via semi-annihilation has been explored in [21].
One of the most crucial parameter, the B −L charge of the scalar dark matter, i.e., q DM is not guided by the model, rather is a free parameter. The non-observation of any direct detection signal motivates to choose a very small value of the charge q DM . Furthermore, a choice of very small q DM suppresses interactions with other particles in equilibrium, leaving out the freeze-out framework completely. In this context some studies have been pursued [9,27,29,30], where B −L gauge boson still contributes in the freeze-in production of the dark matter. The fermion plus scalar dark matter freeze-in scenarios are also explored [31][32][33].
In this article we explore the regime, where the B − L gauge boson contribution is negligible in dark matter relic density, and the freeze-in dynamics is governed by annihilation and decays of SM and B − L scalars. To evaluate the relic density, we adopt the relativistic framework [34][35][36], where we use Bose-Einstien (BE), and Fermi-Dirac (FD) statistics. The effect of thermal mass correction of SM Higgs boson along with electroweak symmetry breaking (EWSB) on dark matter phenomenology has been explored in [37,38] for singlet scalar extension. We explore such effects of EWSB, thermal mass correction, and quantum statistics for scalar extended gauged B − L model. A number of SM and BSM decay and annihilation modes become open at different epoch of the early Universe, which we carefully include in our numerical computation. Similar to [35], which found large enhancement in fusion process, we find significant enhancement in 2 → 2 annihilation, and 1 → 2 decay processes during EWSB, once thermal mass correction of SM Higgs boson has been taken into account. Though our study is confined to gauge B − L framework, the results of this analysis are more generic for mainly two reasons: i) the freeze-in here is dominant by the scalar, and not so by the B − L gauge boson, ii) the relativistic effects at EWSB that we observe are applicable to more generic scenario. The gauged B − L model has also been explored for collider phenomenologies. In the B − L scenario, the RH neutrinos are charged under B − L gauge group and thus they can be produced at the colliders via Z BL unlike the Type I seesaw case [22,[39][40][41][42].
As is evident from the above Lagrangian, the model contains quartic interactions involving dark matter-Higgs, as well as dark matter-S fields, that have major impact in determining the dark matter relic abundance. Other than these particles, the model also contains B−L gauge boson Z BL . See [11][12][13] for detail descriptions of the model. Below, we present a brief dicussion on neutrino masses,  • Gauge boson mass. The additional gauge boson from U(1) B−L is represented by Z BL where the mass of Z BL is generated due to spontaneous breaking of the B − L gauge symmetry, and is given by, In the above g BL represents B − L gauge coupling and the vev of S is denoted by v BL . The √ s = 13 TeV LHC search for a massive resonance decaying into dilepton final states puts a strong lower bound on the Z B−L gauge boson mass, i.e., m Z B−L > 5.15 TeV [43]. For our calculation, we consider m Z B−L = 5.5 TeV, which is in agreement with the LHC bound.
• Scalar masses. Owing to the non-zero λ Sh term in eq. (2.1), and non-zero vev's v, v BL , the scalar fields S and Φ mix with each other after electroweak symmetry breaking. We define the neutral components of S and Φ fields as S + iS I and h + ih I , respectively, which leads to the mass matrix of h and S after EWSB as, Rotating the basis h and S to new states h 1 and h 2 by suitable angle α, we can diagonalise the above mass matrix. The physical mass basis are given by, where h 1 is the SM-like Higgs boson and h 2 is the BSM scalar. The mixing angle between them is given by, (2.5) The mass square eigenvalues of scalar field h 1 and h 2 are given by, For our analysis we stay in the decoupling limit i.e., α ∼ 10 −4 − 10 −5 obeying 2σ constraints of Higgs data of LHC at 13 TeV [44,45]. Therefore, for all practical purposes, due to the very tiny mixing between the SM Higgs and B − L Higgs bosons h 1 h and h 2 S in our analysis. In the subsequent sections, we explore the JHEP05(2021)150 production of dark matter from the SM and B−L Higgs boson decay and annihilation processes. For the above mentioned values of the Higgs mixing angle α and quartic coupling λ Sh > 6 × 10 −6 [35], the B − L Higgs boson is in thermal equilibrium along with SM Higgs boson in the early Universe.
Note that eq. (2.6) represents the physical masses of the scalar fields without any thermal correction and we will see that the thermal correction to the SM Higgs mass have a large impact on the dark matter phenomenology. The electroweak phase transition (EWSB) can be either second order phase transition or a cross-over. The SM Higgs becomes massless in second order phase transition, whereas it remains massive in cross-over during EWSB [37,46,47]. The authors have performed numerical lattice Monte Carlo simulations to study the thermodynamics of the cross-over where they have shown that m h (T c ) approaches around 10-15 GeV during EWSB [47]. In our work, we consider the electroweak phase transition to be a crossover in which Higgs remains massive at critical temperature (T c = 160 GeV). For the temperature is greater than the critical temperature i.e., T > T c , the mass of Higgs bosons is given by [35], whereas for the temperature smaller than the critical temperature i.e., T < T c , the mass of Higgs boson is given by, In the above, c represents a constant which is determined by requiring m h (0) = 125.5 GeV, i.e., Higgs boson mass at zero temperature. In figure 1 we show thermal corrections to the SM Higgs boson mass for different scenarios which we detail later.

JHEP05(2021)150
The vertical lines represent the z = m φ DM T values corresponding to EWSB and for our analysis m h (T c ) ≈ 10 GeV. Note that, due to the difference in the DM mass, EWSB (T EW = 160 GeV) corresponds to different values of z for these different scenarios. This is clearly evident from figure 1. EWSB has a significant importance in our work, as it will be clear from the discussions of the subsequent sections.
Similarly, the thermal correction for the mass of B −L scalar S can also be calculated. However, for our analysis this is not so important and it can be understood easily as follows. To evaluate thermal correction, the parameter µ S in the Lagrangian eq. (2.1) should be replaced by [35], The critical temperature T v BL c is the temperature, where B−L scalar S takes vaccuum expectation value v BL and breaks the U(1) B−L symmetry. The critical temperature T v BL c then can be approximated as [35], In our present work, we consider that the B − L breaking took place at a high temperature T v BL c in the early Universe, which we assume to be equal to the reheating temperature of the Universe [48,49]. Using eq. (2.11) for our parameter choices, we obtain that the re-heating temperature to be T R ∼ 2.26 × 10 4 GeV. Immediately after the re-heating or B − L symmetry breaking, the field S acquires a mass 200 GeV, that we consider throughout our analysis. Hence, thermal correction to S mass is not relevant in our study.
• Neutrino masses. The masses of the light SM neutrinos are generated via the usual Type-I seesaw mechanism: 12) where m N,k = λ N S S are the Majorana masses of the RH neutrinos generated due to the spontaneous symmetry breaking of B − L gauge symmetry.
• Dark matter mass. The mass square eigenvalue of dark matter field φ D is given by, In this work, we consider the couplings λ SD , λ Dh to be very small ∼ 10 −10 -10 −13 to accommodate φ D as non-thermal dark matter. We also consider λ D of similar order 10 −10 , which suppresses any large contribution from 2 → 4 processes, that could have brought the dark matter into kinetic and chemical equilibrium [36]. Due to the JHEP05(2021)150 choice of a small λ D , its impact on the thermal correction of dark matter mass would be negligibly small. This also implies negligible impact of the phase transitions for our choices of dark matter masses which are in the range of a few GeV. To a good approximation, we therefore identify that the dark matter mass is primarily governed by the bare mass term, i.e., m φ DM ∼ µ D and ignore the thermal mass correction of the dark matter.
Before finishing this section, we present a brief discussion about the stability of the dark matter in this model. This is to note, that the dark matter does not acquire a vev in this model. However, since dark matter is charged under B − L and the same symmetry is broken due to non-zero vev of S field, hence the dark matter will not be a stable dark matter for all values of q DM . As given in table 1, the dark matter candidate φ D has charge q DM under U(1) B−L . By choosing appropriate q DM with a value q DM = ±2n (n ∈ Z and n ≤ 4), one can avoid Yukawa interaction terms, such as, φ DN c N and cubic and quartic interaction term, such as φ D S 2 and φ D S 3 [21]. Therefore, the decay of φ D can be forbidden without invoking extra discrete symmetry in the model and hence φ D can be the viable stable dark matter candidate. For dark matter in a different representation other than being SU(2) L singlet, additional re-normalizable and non-renormalizable operators involving SM fields may present, which can further contribute to dark matter decay. This has been studied in [50].
In this work, we consider that φ D is a dark matter with feeble interaction strengths (FIMP candidate). Therefore, in the early Universe, the state had negligible abundance and during reheating of the Universe it was not in the thermal equilibrium. The dark matter φ D has both U(1) B−L gauge and scalar interactions. The production of φ D through gauge interactions are determined by gauge coupling g BL along with the charge q DM of φ D state, the dark matter mass m φ DM and B −L gauge boson mass m Z BL . Here we primarily focus on the dark matter production via the scalar states and for this purpose the q DM is chosen to be sufficiently small, such that, the production of φ D through gauge interactions becomes negligible. In the next section, we present a relative comparison between these two different production modes to justify our choice of parameters.

Freeze-in production of dark matter
As outlined in the previous section, the dark matter particle φ D has feeble interactions with the SM particles, as well as, other B − L particles (S, Z BL ) present in this model. Therefore, the state φ D is not in thermal equilibrium, rather produced from the decays and annihilation of SM and B − L particles. If kinematically allowed, the freeze-in production of dark matter is dominated by the decays of SM and B − L states which are in thermal equilibrium. The production processes due to annihilation give subdominant contributions to the relic density, as often the contributions are suppressed by additional couplings as well as propagators, along with the numerical factors arising from additional phase space integral. A non re-normalizable interactions between the dark matter and bath particles leads to UV freeze-in of dark matter which depends on the re-heating temperature of the JHEP05(2021)150 Universe [1,[51][52][53]. In this work, we do not have non-renormalizable interaction between the dark matter and bath particles. Rather, our scenario is similar to IR freeze-in of dark matter, where production of the dark matter dominates at T ≈ M of the initial states and it is insensitive to the reheating temperature of the Universe. We consider both the decay and 2 → 2 annihilation contributions in the relic density. Feynman diagrams for the production processes of dark matter φ D before and after EWSB are shown in figure 3. Depending on the primary production mechanism, we sub-divide the entire discussion in different Scenarios, and analyse the production of φ D in detail. The schematic diagrams for these different scenarios have been shown in figure 2.
• Scenario-1. The dark matter is primarily produced from the annihilation of the SM Higgs boson.
• Scenario-2. The dark matter is produced primarily from the annihilation of the B − L Higgs, with a sub-dominant contribution from the annihilation of the SM Higgs boson.
• Scenario-3. The dark matter production is governed by the annihilation of the SM Higgs boson at an earlier epoch, but later dominated by the decay of the B − L Higgs bosons.
• Scenario-4. The dark matter production is governed by the decays of SM and B − L Higgs bosons.
• Scenario-5. The dark matter is produced mainly from the decay of SM Higgs boson with a sub-dominant contribution from the B − L Higgs boson. In the earlier epoch, the dark matter production is primarly governed by the SM Higgs annihilation.
We explore each of these different scenarios in detail taking into account all the relevant contributions in the Boltzmann equation. However, before focusing on the main study of this paper, we bring the attention of the readers on a comparative study between the B − L gauge boson (Z BL ) contribution and B − L scalar (S) contribution to the dark matter relic density. It is well known, that when the dark matter is gauged it would quickly thermalise due to potentially larger effective gauge coupling and charge associated with it. In our case such a phenomena can happen as dark matter can be copiously produced via Z BL decay and/or by the annihilation mediated by Z BL or via contact interaction. Such process can lead to overproduction of dark matter in the very early Universe, and the only viable option to maintain the correct dark matter relic is freeze-out [21,22]. However, our goal for this article is to investigate the possibility of relativistic freeze-in scenario which compels us to choose a very small value of q DM . The number density of φ D from Z BL decay i.e., Z BL → φ * D φ D and annihilation processes i.e., ff → φ * D φ D are calculated by the following Boltzmann equation, where Γ aa→bb and Γ a→bb represent the reaction rates for annihilation and decays. In comoving volume the above Boltzmann equation can be written as, The relic abundance of φ D here is mostly dominated by Z BL decay and is given by, Figure 4(a) represents the relic density contours where, we vary the charge q DM and the coupling g BL . This is to note, that the production of φ D via Z BL mediated annihilation processes, i.e.,f f,N c N → Z BL → φ * D φ D are also kinematically allowed but such processes are suppressed due to fourth power of g BL , as well as large Z BL mass compared to Z BL → φ * D φ D process. Therefore, in figure 4(a), we ignore the contributions from the annihilation processes mediated by Z BL . In the same plot, we also show the contour that satisfies the present relic density Ωh 2 = 0.12 [54] by the red dashed line. We can easily infer that for q DM 10 −10 , Z BL → φ * D φ D contribution alone attains the desired relic density. On the other hand, the star mark in figure 4(a) represents our chosen benchmark point, for which the Z BL → φ * D φ D decay gives negligible contribution in the relic density. For a fixed g BL , as we increase q DM , contribution from Z BL decay in relic density will increase. For a very large q DM ∼ 10 −1 , the DM will thermalise with the SM particles, and freeze-out scenario will be the viable option [21].
where, we show the ratio of the relic densities. For our choice of masses, as given in the caption of figure 4, the decay of both B − L scalar S and Z BL into two φ D state are kinematically allowed. The ratio increases significantly with the decrease in g BL and q DM , as can be explained from the following equation,

Masses in GeV Couplings
200 300 250 10 −7 5.0 × 10 −12 6 × 10 −6 0.053 1.6 × 10 −11 We choose q DM ≈ 10 −12 represented by the red star in figure 4(a) where it is evident that the production of φ D through gauge interaction is negligible and thus we neglect contribution from the B − L gauge interaction in the rest of the paper. Even if we consider Z B−L mass different from 5.5 TeV [43], as long as we are choosing a sufficiently small q DM , production of dark matter from the scalar sector will continue to dominate. We focus on the production of the dark matter from decay and annihilation of the scalars in relativistic freeze-in scenarios. The effect of SM fields (fermions, gauge bosons) are also taken into account via the interactions which are operative after EWSB. As discussed in the previous section, for our analysis we consider that the re-heating temperature of the Universe is same as the temperature at which B − L symmetry breaks down. To evaluate dark matter number density, we therefore perform the analysis in the B − L broken phase. Below, we present a detailed discussion of the different scenarios, where we numerically solve the Boltzmann equation and evaluate the relic density. In doing so, we consider different decay a → φ * D φ D and annihilation/co-annihilation processes aa, ab → φ * D φ D .

Scenario-1
In this scenario, the dark matter production primarily occurs via SM Higgs boson annihilation i.e., hh → φ * D φ D . We adopt a relativistic freeze-in framework for the evaluation of the relic density. The contribution of SS → φ * D φ D is although allowed but small in attaining the correct dark matter relic. Since m φ DM > m h /2, m S /2 the decay contributions from the SM and B − L Higgs bosons are absent. The choices of masses and coupling used in the numerical analysis, are shown in table 2. The Boltzmann equation for the production of φ D in this scenario is given by eq. (3.5), where Γ aa→bb and Γ a→bb are the rates of the annihilation and decay processes for the respective channels. In a comoving volume the above Boltzmann equation can be read in terms of the yield as described in eq. (3.6),   Figure 5. The figures correspond to Scenario-1. Figure 5(a) shows the relativistic reaction rates for the process hh → φ † D φ D and other relevant processes. Figure 5(b) shows the individual contributions to the relic density, and the total relic density. The brown horizontal line represents the present experimentally measured relic density [54]. where the number density n is related with the entropy density s as n = Y /s and H is the Hubble's constant. Before EWSB, all four degrees of freedom of the SM Higgs doublet contribute to ΦΦ → φ * D φ D leading to four times enhancement in the relic density as compared to hh → φ * D φ D after EWSB (green line). The expressions of the different reaction rates are given in appendix A, where we have closely followed the approach of [34,36]. The annihilation processes hh, SS, N N are always open while the other SM annihilation processes unlatch only after EWSB.
In figure 5(a) we show the relativistic rates for different annihilation processes hh, SS,
• The Higgs annihilation processes ΦΦ, hh → φ * D φ D , 1 is mostly dominated by the contact four point diagram, given in figure 2. However, for our numerical analysis contributions from all relevant diagrams (mediated via h, S) are taken into account. The cross-sections for the processes are listed in appendix B. It is worth mentioning that due to the choice of the dark matter mass none of the annihilation process contains any resonant production. The t-channel diagram gives negligible contribution, and hence, is not considered.
• Similar to the previous case, for the annihilation channel SS → φ * D φ D , dominant contribution arises from the contact term.
• The annihilation rate of N N → φ * D φ D is much suppressed as compared to ΦΦ, hh → φ * D φ D due to additional couplings with S and the corresponding propagators.  5(c) depicts the relative enhancement in the relic density obtained using BE distribution as compared to MB distribution, which is as significant as ∼ 1.42-1.62. The relative interaction strength Γ BE Γ MB also varies accordingly with z.

Scenario-2
Unlike the previous scenario, the φ D production is governed primarily by annihilation of The larger production from B − L scalar annhilation occurs due to a larger λ SD compared to λ Dh , as can be seen from table 3. Similar to the Scenario-1, here also h and S decays are kinematically forbidden. The Boltzmann equation in this case would be the same as eq. (3.5) and so is the yield equation i.e., eq. (3.6). Similar to the previous case here also all the contributions coming from hh, SS, N N → φ * D φ D are taken into account for the numerical analysis. The results are manifested in figure 6. In figure 6(a), we show 1 Φ is the SM Higgs doublet and h represents the SM Higgs field after EWSB. In figure 5(a), we do not maintain this distinction, rather represent the SM Higgs doublet(before EWSB) and SM Higgs field(after EWSB) by h only. We show the contribution from one massive degree of freedom in hh → φ * D φD. In figure 5(b), all the four contributions (before EWSB) and one contribution (after EWSB) have been considered.  Figure 6(a) shows the relativistic reaction rates for different process corresponding to production of φ D . Figure 6(b) shows the individual contributions to the relic density, and the total relic density. The brown horizontal line represents the present experimentally measured relic density [54]. Figure 6(c) shows the relative enhancement in the relic density with respect to Maxwell-Boltzmann distribution.

JHEP05(2021)150
the relativistic rates for different processes corresponding to production of φ D . It can be seen that SS → φ * D φ D (orange line) is the most-dominant, where this channel is governed by the contact interaction. hh → φ * D φ D rate (green line) is the second dominant and its contribution is only 16% to the dark matter relic. Due to the choice of the mass of dark matter, neither of the above two channels entail any resonance enhancement.
In figure 6(b), we show the production of φ D from different processes. In the present scenario, the dark matter freeze-in occurs at the temperature of 150 GeV. We can see that the N N → φ * D φ D annihilation mode is suppressed due to additional couplings with S, Z BL and the corresponding propagators. The SM annihilation to dark matter, i.e., W + W − , ZZ, tt, bb → φ * D φ D starts only after EWSB and are mediated via off-shell h, S. These processes contribute ∼ 3.3% only. Figure 6(c) depicts the relative enhancement of relic abundance using BE distribution over the MB distribution. We can see that the dark matter production is dominated by annihilation and the relative enhancement in the relic density using BE distribution is quite significant ∼ 2.5-1.65.

Scenario-3
Along The decay h → φ * D φ D however still remains forbidden. The Boltzmann equation contains the decay contribution as well, and can be written as, The yield equation in the comoving volume is given by, As before, here also we consider all possible annihilation processes hh, SS, N N , hS → φ * D φ D and decay S → φ * D φ D , along with other SM processes. The results are summerised in figure 7. Figure 7(a) represents the relativistic rates for the different processes corresponding to the production of φ D . Figure 7(b) shows the evolution of φ D . The annihilation mode hh → φ * D φ D (green line) remains most efficient in the production of the dark matter for z < 0.02. After this S → φ * D φ D (purple line) takes over and remains the dominant mode until dark matter freezes-in at the temperature of 40 GeV. Since the S decay is open, an abrupt increase in λ SD can cause an overproduction of φ D . Around EWSB where the Higgs mass falls below 100 GeV, s-channel resonance occurs in the process hh → φ * D φ D (green bump) which enhances the production rate.
Here we elaborate on an important point about our calculation which resolves the possible over counting. We take into account all possible contributions in the evaluation of hh → φ * D φ D process, namely the contact term, on-shell and off-shell contributions mediated by h, S. The S mediated s-channel diagram encounters a resonance around EWSB. In this case the production of φ D must be computed by subtracting the on-shell S exchange contribution to avoid any over-counting, as in the Boltzmann equation this effect has already been taken into account by the S decay contribution separately (eq. (3.7)) [55].  Figure 7(a) shows the relativistic reaction rates for different process corresponding to production of φ D , and S → φ D φ D . Figure 7(b) shows the individual contributions to the relic density, and the total relic density. The brown horizontal line represents the present experimentally measured relic density [54]. Figure 7(c) shows the relative enhancement in the relic density with respect to Maxwell-Boltzmann distribution.

JHEP05(2021)150
The on-shell contribution due to S mediation is given by, Therefore, only contact term and off-shell contribution hh are taken into account for the φ D production in this scenario. Figure 7(c) depicts the relative enhancement of relic abundance using BE and MB distribution. One can see that at a very early epoch, where the dark matter production was dominated only by SM Higgs boson h annihilation, the ratio is very significant around ∼ 1.6. At the later epoch, when the production of φ D is dominated by the S decay, we find that the enhancement is about 1.04.

Scenario-4
This is the most generic scenario where along with different annihilation processes, the decays of both the Higgs bosons h, S → φ * D φ D are kinematically allowed. The chosen

Masses in GeV Couplings
200 300 1 10 −7 6.65 × 10 −13 6 × 10 −6 0.053 8.6 × 10 −12 benchmark points are tabulated in table 5. Unlike previous cases, the dark matter in this scenario is very light m φ DM = 1 GeV. As we will show in the subsequent discussion, the decay of Higgs bosons h, S give the most dominant contribution in the relic density. The most generic Boltzmann equation involved in this case has the following form: (3.10) The corresponding yield evolution in the co-moving volume can be written as, We include all the contributions as mentioned in eq. (3.10) in our numerical analysis. We show the relativistic reaction rates, relic density and the relative enhancement of relic density in figure 8. In figure 8(a), we illustrate the relativistic reaction rates for different processes corresponding to the production of φ D . Similar to Scenario-3, here also we observe the resonant enhancement around EWSB for hh → φ * D φ D annihilation (green bump). The reaction rate for S → φ * D φ D dominates until almost EWSB, after which h → φ * D φ D takes over. Similar to the previous scenario, we avoid any over-counting of S on-shell production, by removing it from hh → φ * D φ D annihilation process. Such procedure has also been followed for other similar processes, such as, bb  Figure 8(a) shows the relativistic reaction rates for the process hh, SS → φ D φ D , and h, S → φ * D φ D . Figure 8(b) shows the individual contributions to the relic density, and the total relic density. The brown horizontal line represents the present experimentally measured relic density [54]. Figure 8(c) shows the relative enhancement in the relic density with respect to Maxwell-Boltzmann distribution.

Scenario-5
In this scenario, the primary contribution to relic density arises from the decay process h → φ * D φ D as we choose higher λ Dh . Similar to Scenario-4, here also we choose a light dark matter with mass m φ DM = 1 GeV, as shown in table 6. The Boltzmann and the yield equations have the same form as in Scenario-4, so we follow eq. (3.10) and eq. (3.11) for our numerical analysis. Figure 9 shows variation of the different annihilation and decay  Figure 9. The figures correspond to Scenario-5. Figure 9(a) shows the relativistic reaction rates for the process hh, SS → φ * D φ D , and h, S → φ D φ D . Figure 9(b) shows the individual contributions to the relic density, and the total relic density. The brown horizontal line represents the present experimentally measured relic density [54]. Figure 9(c) shows the relative enhancement in the relic density as compared to Maxwell-Boltzmann distribution and figure 9(d) shows the relative enhancement in the respective reaction rates. channels along with the evolution of dark matter relic and BE/MB comparison wih z. Figure 9(a) shows that the annihilation hh → φ * D φ D is dominant only at a very early epoch, then the S decay takes over, and finally at EWSB i.e., z ≈ 6.25 × 10 −3 , h decay opens up and becomes the most dominate till φ DM freezes-in around T ∼ 20 GeV. Due to the choice of a light dark matter, the process hh → φ * D φ D mediated by S encounters schannel resonance during EWSB, which is shown by the green bump. We follow the same prescription as before, where we omit on-shell contribution from the above mentioned process, and consider only contact term and off-shell contributions in the relic density. Other SM annihilation processes, such as, bb, tt, W + W − , ZZ → φ * D φ D (mediated via h, S) open up only after EWSB, however, their rates are relatively small. Similar to the previous scenario, the bb → φ * D φ D also entails a resonance, due to mediation of an on-shell h, and we again adopt the same prescription as Scenario-4. Since the N N contribution is very small, we do not show that in figure 9(a), and figure 9(b). None of the any other SM annihilation channels contain any resonance.

JHEP05(2021)150
In figure 9(b) different contributions in obtaining the correct dark matter relic are shown. It is seen that h → φ * D φ D and S → φ * D φ D are two dominant modes with 83% and 17% contributions towards attaining the desired dark matter relic in a freeze-in mechanism at the temperature of 20 GeV. h → φ * D φ D is the leading contributor due to a larger λ Dh as compared to λ SD . The other SM contributions bb, tt, W + W − , ZZ → φ * D φ D are small ∼ 0.1% only. Figure 9(c) depicts the relative enhancement of relic abundance using BE and MB distributions as a function of z. At an early epoch, hh → φ * D φ D primarily dominates the relic density leading to a ratio 1.6. The ratio falls as temperature decreases and at EWSB, i.e., z ≈ 6 × 10 −3 a distinct kink appears in the ratio. The kink is more pronounced as compared to Scenario-4, due to the presence of SM Higgs decay h → φ * D φ D at EWSB. This kink is due to the sudden jump in the rates as can be seen from figure 9(d). This is also to note that for the annihilation hh → φ * D φ D (via contact) (green line) kink is more pronounced than the decay process h → φ * D φ D (red line) at EWSB. The relative enhancement in the relic density varies from ∼ 1.6 at a lower value of z to ∼ 1.02 during the freeze-in temperature of 16.66 GeV.
The nature of the ratio of the relic densities obtained from BE and MB distributions and the appearance of the kink can be understood in the following way. The 1 → 2 rate has the following form: For SM Higgs boson (h) decay using BE distribution we get (3.14) Whereas for MB distribution this rather becomes, Now we can compare them as follows .

(3.16)
At EWSB when T = 160 GeV the Higgs mass becomes m h = 10 GeV (see figure 1). At that point the next to leading order terms contribute substantially to the Γ BE However, this is not true for S decay, since mass of S is considerably large m S = 200 GeV, which does not lead to this kind of enhancement at EWSB.

JHEP05(2021)150
Overall, we find that the thermal mass correction to SM Higgs boson and the freeze-in temperature have a large impact in determining the final enhancement factor. For Scenario-1,2 the freeze-in occurs around EWSB. The hh → φ * D φ D rate is substantially large around EWSB (this holds for all scenarios), when BE distribution is being used. This leads to a large final enhancement in the ratio of relic density Ωh 2 BE /Ωh 2 MB . However, for Scenario-3,4,5 the freeze-in occurs at a later epoch than EWSB and hence the final enhancement factor is relatively small. This phenomena is rather generic and one can see it from the first term of eq. (3.14) which corresponds to reaction rate using MB distribution (see eq. (3.16)) and other terms in the series provide the correction to the reaction rate. The behaviour of K 1 (z) plays significant role in determining the relative enhancement in the reaction rates and hence the relative enhancement in the relic densities. Whenever decaying particle's mass becomes less than temperature i.e., m h T , z 1, the correction gives significant contribution to the rates. This occurs in the limit where the Bessel function i.e., K 1 (z) ≈ 1/z 1. Similarly, for m h T i.e., z 1 limit, the Bessel function K 1 (z) ≈ e −z √ z 1 [56], which results into the same reaction rate for both BE and MB distributions.

Conclusion
We analyse the freeze-in production of a scalar dark matter in an extended gauged B − L model where a complex scalar field φ D is the dark matter candidate. To evaluate its relic abundance, we follow a relativistic formalism, where we consider Bose-Einstein and Fermi-Dirac statistics. Due to a very tiny charge of the dark matter under B −L gauge symmetry, that in turn leads to a suppressed interaction of dark matter with the B − L gauge boson, its production from the Z BL gauge boson is negligible, and hence not important for our study. We rather focus on the annihilation and decay of the SM and B − L Higgs boson hh, SS → φ * D φ D , and h, S → φ * D φ D that contribute primarily to the relic density. In evaluating the annihilation contribution, we consider all possible processes, namely, the contribution from the four-point contact interaction involving Higgs/B − L Higgs boson and dark matter, that directly contribute to hh/SS → φ * D φ D , as well as, any other schannel mediated processes. The t-channel diagrams give suppressed contribution in our case, and hence have not been considered. Additionally, we also consider the annihilation of SM particles, such as W + W − , ZZ, bb, tt → φ * D φ D that contribute at most by 1% to the relic density. Depending on the mass of dark matter, and the primary production mechanism, we consider five different scenarios Scenario 1-5. We show that thermal correction to the SM Higgs mass has a significant impact on different annihilation channels, such as, where, a few of these processes undergo resonance enhancement in their respective reaction rates, due to the on-shell mediation of h, S states. We consider a fixed mass m S = 200 GeV for this study. The entire discussion have been sub-divided into the following few scenarios.
• In Scenario-1 and Scenario-2, we explore the freeze-in production assuming a dark matter with mass 250, 150 GeV, respectively. The primary dark matter production mechanism is the hh → φ * D φ D for the 1st scenario, and SS → φ * D φ D for the second JHEP05(2021)150 scenario. Due to the choice of the mass of dark matter, neither the SM or B −L Higgs boson decays to dark matter state. We find for a large λ Dh coupling, the hh → φ * D φ D dominates the production, whereas for a large λ SD , the SS → φ * D φ D gives dominant contribution.
• In Scenario-3, we consider dark matter mass to be 80 GeV. This serves as one of the illustrative cases, where both the annihilation of SM Higgs boson hh → φ * D φ D , and the decay of B − L Higgs boson S → φ * D φ D contribute to the relic density. The SM Higgs annihilation channel serves as a primary production channel at an early epoch, while the S → φ * D φ D channel becomes dominant at a later epoch. Due to choice of dark matter mass, the Higgs boson decay is kinematically forbidden.
• In Scenario-4 and Scenario-5, we consider the dark matter mass to be significantly lower than the SM Higgs boson mass, m φ DM = 1 GeV. For Scenario-4, both the h, S decays contribute almost equally to the relic density. For Scenario-5, the primary production mode is the SM Higgs boson decay to dark matter particle.
This article presents a comparison between the relic density obtained by using BE statistics, with the one obtained by using MB statistics. We see for the annihilation dominated scenarios, Scenario-1,2, where freeze-in occurs during EWSB, the final ratio of relic density obtained using BE and MB statistics is large, R = Ω BE h 2 Ω MB h 2 varies between 1.42-1.62. For the other three scenarios, where the decay of SM and B − L Higgs bosons dominate the relic density and freeze-in occurs at a much later epoch for the Scenario-3,4,5, the enhancement factor is much less 1.04.
This effect is inherently linked with thermal mass correction of SM Higgs boson, which is considered in this study. However for the scenarios considered here, the thermal mass correction to S and dark matter are not relevant. We consider the EWSB as a crossover and explore the effect of thermal mass correction of SM Higgs boson on dark matter abundance. It is noticed that due to the low mass of the SM Higgs boson i.e., m h = 10 GeV at EWSB temperature T = 160 GeV, the relativistic reaction rate for hh → φ * D φ D via contact term becomes significantly enhanced. This occurs as the correction terms in the relativistic reaction rate obtained using BE statistics become significantly large during EWSB. This results in an enhanced reaction rate of hh → φ * D φ D during EWSB, when BE statistics being used in particular relevant for the light dark matter mass. The relative enhancement is more pronounced for annihilation (almost 2.3), as compared to decay ( 1.5).
The relative enhancement in the reaction rates also result in a distinct kink in the R around EWSB for Scenario-4,5, where SM Higgs boson decay or annihilation processes contribute significantly in the dark matter relic abundance. For Scenario-2,3 since the S decay or annihilation are dominant, therefore, we do not see such an intermediate kink in the ratio of relic density plot.
We conclude with the observations that quantum statistics, along with the thermal mass correction are essential to capture these enhancement effects in dark matter relic density in freeze-in scenario which otherwise would be overlooked.

JHEP05(2021)150
Finally, we make qualitative remarks about few other possible extensions. Allowing a large dark matter self-interaction, its coupling with the SM and BSM Higgs field, and also a much suppressed gauge coupling g BL will lead to a different freeze-in dynamics. In this work, we have assumed dark matter self-interaction is negligible. For a large dark matter self-interaction, dark matter will thermalise with itself via 2 → 4 processes in the early Universe [36,57,58]. We have also considered, that the quartic interactions of the dark matter with SM and BSM Higgs is negligible. For our assumptions about the dark matter self-interaction, and quartic coupling with the scalar fields, dark matter in our scenario is non-thermal. Furthermore, due to small couplings associated with quartic interactions their impact on the thermal correction of dark matter mass is negligible. Allowing a large quartic coupling, the thermal correction to the dark matter mass will be sizeable, which needs to be included in the study. For a large dark matter self-interaction, 2 → 4 processes will also be important.
Furthermore, for a suppressed g BL coupling, freeze-in dynamics will be much more involved. For our chosen benchmark points, which includes a large gauge coupling g BL , and the charge of BSM Higgs, the BSM Higgs S quickly thermalises with the SM particles. Hence, sequential freeze-in [59], i.e., production of S and then production of the dark matter from S can not be materialised. However, for a very suppressed gauge coupling g BL by many orders of magnitude, and also suppressed quartic interactions of the BSM Higgs with SM Higgs field, most of the BSM (BSM Higgs, heavy neutrino N , Z BL etc) particles in our model will be non-thermal. Their evolutions in the early Universe will be highly dynamic and coupled which will be determined by solving sets of coupled Boltzmann equations. The detail analysis of these few interesting possibilities is beyond the scope of this paper, and will be explored in a further study.

JHEP05(2021)150
been derived in [34,36]. Here, we briefly summarise the results. The reaction rate per unit volume has the generic expression: (A.1) Here M a→b is the transition amplitude and f (p) is the momentum distribution function.
In thermal equllibrium, f (p) can be written in a covariant form as where the upper(lower) sign is for fermionic (bosonic) particles. The final states can either be in equilibrium or non-equilibrium with thermal bath. For the final states, not in equilibrium with the thermal bath implies a negligible initial abundance, leading to the final state enhancement factor 1 + f (p j ) ≈ 1. Similarly, for the final states which are in equilibrium with thermal bath, one can neglect the pauli-blocking /stimulated emission effects, i.e., 1 + f (p j ) ≈ 1. For 2 → 2 processes, cross section is defined by The reaction rate can be written in terms of cross section which is given by where v mol is the moller velocity of the incoming particle, and is given by, The reaction rate can be easily evaluated in centre of mass (CM) frame. See [34,36] for the details. Following [34,36], we define two new variables p = (p 1 + p 2 )/2 and k = (p 1 − p 2 )/2 for a pair of momenta p 1 and p 2 . The vector p can be Lorentz transformed to the form In the above, E represents the particle energy in CM frame. In terms of half of the centre of mass energy E, rapidity η and angular coordinates θ,φ, the vector p can be expressed as [34] p 0 = E cosh η,

A.1 Annihilation
The reaction rate for ab → cd processes of incoming bosons is given by . We derive the reaction rates for 2 → 2 processes of incoming fermions which has the following expression: At the low temperature T limit, Therefore, at the low temperature T limit, eq. (A.8) and eq. (A.9) reduces to the reaction rate which is equal to the reaction rate obtained by MB distribution, For the incoming particles having same mass m = m a = m b , eq. (A.8) and eq. (A.9) reduce to and (A.12) for the incoming fermions. At low temperature T limit, the eq. (A.11) and eq. (A.12) reduces to the reaction rate which is equal to the reaction rate obtained by MB distribution, At the low temperature T , eq. (A.20) and eq. (A.21) reduces to a simpler form which is equal to the decay rate obtained by MB distribution, Next to the leading order terms contribute in eq. (A.20) and eq. (A.21) substantially when M T . Therefore, the Γ BE 1→2 is enchanced and Γ F D 1→2 is suppressed when compared with Γ MB 1→2 .

B Analytical expressions of relevant cross sections and decay widths
We provide the expressions for the relevant cross sections in half of the centre of mass frame and decay widths, that have been used in the Boltzmann equations.

B.1 Cross sections for different processes
Below we consider that the B − L Higgs boson S, and SM Higgs boson h, and we neglect the mixings between SM and B − L Higgs boson. As considered in the text, the SM and B − L Higgs boson mixing angle sin α = 10 −4 − 10 −5 .
The first term is from contact term λ DH Φ † Φφ † D φ D . The second and third terms are via mediation of S and h respectively, where third term appears only after EWSB.

SS
The first term is from contact term λ DS S † Sφ † D φ D . The second and third terms are via mediation of h and S respectively, where second term appears only after EWSB.
We consider the S mediated diagram, and ignore the Z BL , h mediated processes. For Z BL , the Z BL − φ D − φ D coupling is vanishingly small. Since we consider the SM and B − L Higgs boson mixing to be very small, hence, we ignore the h mediated diagram.

JHEP05(2021)150
This process is mediated by h, S. The W W, ZZ, ff → φ † D φ D processes are mediated via only h. Since the SM and B − L Higgs boson mixing angle is tiny, we neglect the S mediated contribution.
In the above, n c is the color charge, and is 1 for leptons and 3 for quarks.

B.2 Decay widths of S
The expressions for the decay widths of S: In the above, n c represents the color charge, and is 1 for leptons, and 3 for quarks.

C Thermal correction to SM Higgs mass
In this work, we have considered the electroweak phase transition to be crossover in which Higgs remains massive at critical temperature (T c = 160 GeV). We have assumed m h (T c ) ≈ 10 GeV. When the temperature is greater than the critical temperature T c . In this regime, the mass of Higgs bosons is given by [35], where c is a constant determined by the requirement m h (0) = 125.5 GeV.
Open Access. This article is distributed under the terms of the Creative Commons Attribution License (CC-BY 4.0), which permits any use, distribution and reproduction in any medium, provided the original author(s) and source are credited.