Secluded Dark Matter in Gauged $B-L$ Model

We consider the gauged $B-L$ model which is extended with a secluded dark sector, comprising of two dark sector particles. In this framework the lightest $\mathcal{Z}_2$-odd particle is the dark matter candidate, having a feeble interaction with all other SM and BSM states. The next-to-lightest $\mathcal{Z}_2$-odd particle in the dark sector is a super-wimp, with large interaction strength with the SM and BSM states. We analyse all the relevant production processes that contribute to the dark matter relic abundance, and broadly classify them in two different scenarios, a) dark matter is primarily produced via the non-thermal production process, b) dark matter is produced mostly from the late decay of the next-to-lightest $\mathcal{Z}_2$-odd particle. We discuss the dependency of the relic abundance of the dark matter on various model parameters. Furthermore, we also analyse the discovery prospect of the BSM Higgs via invisible Higgs decay searches.


Introduction
Observational evidence shows that 84% matter of the universe is in the form of non-baryonic dark matter. However, very little is known about the nature of dark matter (DM) and its origin. The Standard Model (SM) can not explain the observed relic density. One of the most favoured scenarios for DM production has been thermal freeze-out, where a weakly interacting massive particle (WIMP) serves as a DM candidate. The WIMP with electroweak-scale mass, which interacts with other particles via electroweak interaction, can naturally explain the measured DM relic density of Ωh 2 = 0.1199 ± 0.0027 [1]. Several direct detection experiments so far have searched for a WIMP. However, the lack of conclusive experimental evidence motivates the exploration of alternate dark matter production mechanisms. The production of DM via freeze-in mechanism [2][3][4][5][6][7][8][9][10] is one of the most popular production mechanisms. In this framework, DM has a very tiny interaction with the SM and any other particle which are in thermal equilibrium with the plasma and thereby referred to as feebly interacting massive particle (FIMP). Due to significantly suppressed interaction with SM/BSM particles, the FIMP DM never attains thermal equilibrium. Due to a similar suppression in the interaction with the SM particles, FIMP, in general, can not produce any observable signal in the direct detection experiments. See [11] for other alternatives with a light mediator. The DM in the freeze-in scenario is produced from the decay and/or annihilation of SM, and BSM particles which are either in equilibrium with the thermal plasma or also freezing-in along with the DM [12].
Apart from the DM abundance, the SM fails to explain neutrino masses and mixings. One of the most promising models that explain small neutrino masses is the gauged B − L model, which contains three right-handed neutrinos (RHNs) [13][14][15][16] that generate light neutrino masses via seesaw mechanism [17,18]. In addition, the model also contains one BSM gauge boson Z BL and a complex scalar field S. The scalar field acquires vacuum expectation value and breaks the B − L gauge symmetry. The BSM gauge boson and the heavy neutrinos acquire their masses due to the spontaneous breaking of the B − L gauge symmetry. The model can further be extended with a secluded dark sector with a scalar particle φ D with non-zero B − L charge and a gauge singlet fermion state χ, where either or both of them can be suitable DM candidates. The dark sector particles are odd under a Z 2 symmetry. The thermal DM for this model has been explored in several works, such as [19,20]. For a different variation of the B − L gauge model with only a thermal scalar DM, see [21][22][23][24][25]. For RHN DM in a typical B − L model, see [26][27][28][29][30][31]. The RHN can also serve as a portal between the SM sector and a secluded dark sector containing DM particles, see [32][33][34] for the relevant discussion. One of the interesting possibilities is if the fermion state χ serves as the non-thermal DM candidate, and the scalar particle φ D , which was in thermal equilibrium in the early universe, has a significant contribution in the production of χ.
In this article, we consider this possibility, where χ is the DM state, and φ D significantly impacts its production. We study the production of DM through a thermal freeze-in mechanism and significant non-thermal freeze-in contribution [3,35] from the decay of φ D . The state χ interacts only with φ D and RHN N . The DM candidate χ never thermalises due to very tiny coupling. This work considers that φ D is heavier than χ state. However, it serves as the next-to-lightest Z 2 -odd particle (NLOP). φ D thermalises with the SM particles because of large quartic couplings in the scalar potential, sizeable SM-BSM Higgs mixing angle, and large gauge coupling. DM is primarily produced at high temperatures via a thermal freeze-in mechanism from the decay of bath particle φ D , which was in equilibrium with the rest of the plasma. φ D subsequently decoupled from the thermal bath, and its late decay further produced substantial DM relic density. The abundance of φ D at the time of decoupling Ω F O φ D h 2 is governed by the freeze-out mechanism. Depending upon the abundance Ω F O φ D h 2 and its conversion to χ state through out-of-equilibrium decay, the DM can primarily be produced from the late decay as well, which we refer as non-thermal production of DM. We divide the discussion into two different scenarios and show the importance of thermal and non-thermal freeze-in contributions in determining DM relic abundance. We elaborately discuss the constrain on other model parameters, such as the scalar quartic couplings, SM-BSM Higgs mixing angle and the mass of φ D and χ, that appear from the relic density constraint. Other than this, we also explore the discovery prospect of this model at the future High Luminosity run of the LHC (HL-LHC), mainly focusing on the heavy Higgs searches. Due to non-zero coupling λ SD between the state φ D and B − L Higgs S, the model offers a largely invisible branching ratio of the BSM Higgs state. We analyse the discovery potential of the BSM Higgs H 2 -φ D scalar quartic coupling at the HL-LHC in vector-boson fusion (VBF) channel. This quartic coupling has a large impact in determining φ D abundance, hence the DM relic density.
The paper is organised as follows. In Section. 2, we describe the model. Following this in Section. 3, we discuss the dark matter production in detail and along with the effect of out-of-equilibrium decay of φ D . And later, we discuss the collider prospects on the search of φ D at the LHC in Section. 4. We present our conclusion in Section. 5. In Appendix. A.1, A.2 and A.3, we provide the necessary calculation details.

Model
The model is a gauged B − L model, augmented with a secluded dark sector. In addition to the particles of the gauged B − L model, i.e., the right-handed neutrinos N , BSM Higgs field S and BSM gauge boson Z BL , the model also contains dark sector particles a complex scalar state φ D and a B − L singlet fermion χ. The dark sector particles are odd under Z 2 symmetry, while other particles are even. The Z 2 symmetry ensures the stability of DM. We consider that the BSM sector of the gauged B − L model and the SM Higgs state h act as a portal between other SM particles and the dark sector 1 .
The charge assignments of different particles are shown in Table 1. Here, the field S represents a complex scalar field, which acquires vacuum expectation value (vev) v BL = 0 and breaks B − L gauge symmetry. The state N contributes to the light neutrino mass generation via the seesaw mechanism. Note that the scalar φ D is non-trivially charged under B − L gauge symmetry, while the fermion χ is a singlet under both the B − L and SM gauge group. As we will see in the subsequent sections, this leads to significant differences in the evolution of χ and φ D abundances. The complete Lagrangian of the model has the following form, where L DM is the Lagrangian containing dark sector particles, and L B−L is the B − L Lagrangian. The B − L Lagrangian has the following form, 2) 1 In our notation, h represents the SM Higgs doublet field.
The dark sector Lagrangian is given by, The interaction strength between φ D and BSM and SM Higgs bosons(i.e, S and h) are proportional to λ SD and λ Dh . The kinetic energy terms involving S, φ D and N contains the covariant derivatives which is given by, where X = S, N, φ D and Y B−L (X) represents B − L charge of the states shown Table 1.
• SM Higgs and BSM Higgs: After spontaneous symmetry breaking (SSB), the SM Higgs doublet h and BSM scalar S is given by, Owing to the non-zero λ Sh , h 1 and h 2 mixes with each other which leads to the scalar mass matrix given by, The basis states h 1 and h 2 can be rotated by suitable angle θ to the new basis states H 1 and H 2 . The new basis states represents the physical basis states which are given by, where H 1 is the SM like Higgs and H 2 is the BSM Higgs. The mixing angle between the two states is defined by, The mass square eigenvalues of H 1 and H 2 are given by, (2.10) • Neutrino mass: The Majorana mass term of RHN's is generated through spontaneous symmetry breaking of B − L symmetry. The mass of RHN's is given by, The mass of the SM neutrinos is generated through Type-I seesaw mechanism where the light neutrino mass matrix has the following expression, (2.12) • Gauge boson mass: Similar to the RHN's, the additional neutral gauge boson mass Z BL is generated via spontaneous breaking of B − L gauge symmetry. The mass of Z BL is related to the symmetry breaking scale v BL as, where g BL is the associated B − L gauge coupling constant.
• Dark sector constituents mass: The mass square of the particle φ D has the following form, (2.14) Note that, both electroweak symmetry breaking vev v and the B − L symmetry breaking vev v BL have impact in determining the mass of φ D . In this work, we consider λ SD and λ Dh in between 10 −1 and 10 −5 to have φ D as the thermal particle.
• The DM candidate χ is singlet under SM and B −L gauge group. It's mass is governed by the bare mass term, i.e, m χ .

Dark Matter
The dark sector fields χ and φ D can be DM particles. However, we consider the scenario, where χ is a non-thermal FIMP DM, and φ D , the NLOP, is primarily responsible for DM production. In the early universe, the state φ D was in thermal equilibrium with the bath particles, and at some later epoch denoted as T d , it decoupled from the rest of the plasma. Other than the standard thermal freeze-in contribution via φ D → χN process effective up to epoch T ∼ m φ D > T d , the out-of-equilibrium decay of φ D into χ also contributes significantly in the relic abundance of DM. Below, we explore this possibility in detail.

Super-WIMP φ D + FIMP DM χ
This is to note that the particle χ has only one portal interaction Y Dχχ φ D N with the dark sector field φ D and the RHN field N , where Y Dχ is the respective coupling. We consider Y Dχ ∼ O(10 −10 − 10 −12 ) to be very small, and hence, χ has feeble interaction with every other particle of this model. Due to tiny interaction, it fails to thermalise with the rest of the plasma. It was produced from the decay and annihilation of the SM and BSM particles in the early epoch. We show Feynman diagram for all possible decay and annihilation processes that contribute to the χ production in Fig. 1. Among all these processes, the production of χ from decay processes, however, dominates. The sub-dominant contribution to the production of χ from the annihilation of SM and BSM states arises due to additional small couplings, heavy propagators, and suppression factors arising from phase space integral. This work sticks to the renormalizable interaction between the dark sector, the SM, and B − L particles. Therefore, the production of χ is insensitive to the reheating temperature, set by reheating/end of inflationary dynamics. The abundance of χ builds up primarily due to the φ D → χN process and increases when T > m φ D . The production of χ is most significant when T ≈ m φ D . When the temperature falls below m φ D , i.e., T < m φ D , Boltzmann suppression of the parent state φ D occurs and DM χ freezes in. This is referred to as thermal freeze-in production, as the parent particle φ D during this epoch has been in thermal equilibrium. In addition to the standard thermal freeze-in contribution, the abundance of DM χ can further be enhanced from the out-of-equilibrium decay of φ D . This occurs at a later epoch T << m φ D , when φ D is in out-of-equilibrium. The production of DM from the late decay of a state which has decoupled from the thermal plasma is referred to as the super-wimp mechanism, which has been discussed in [3,[35][36][37][38]. In our scenario, the state φ D serves as a super-wimp candidate. The state φ D having non-zero B − L charges and non-zero quartic couplings λ SD and λ Dh interacts abundantly with the SM and B − L particles. At an earlier epoch, φ D hence was in thermal equilibrium with surrounding plasma, maintaining an equilibrium distribution. The non-thermal decay of φ D , which enhances the DM relic density, takes place at a late stage in the thermal history. The dark sector state φ D tracked equilibrium abundance when the temperature of the universe was greater than its mass. The dark sector states φ D abundance decreases mainly through annihilation which are efficient up until m φ D T ≈ 25. Feynman diagram for φ D depletion through annihilation to B/SM particles is shown in Fig. 2. Around this temperature, the interaction rate for the annihilation/scattering of φ D becomes less than the expansion rate of the universe. Hence, φ D fails to scatter with surrounding plasma constituents, and it decouples from the cosmic soup. To compute the relic density of dark sector constituents, one needs to study the evolution of the number density of its constituents with the temperature of the universe. The evolution of the number density of the dark sector constituents is governed by the Boltzmann equations, which contain all the information of the number changing processes of the dark sector constituents. The Boltzmann equations for the evolution of φ D and χ in terms of its co-moving number density Y φ D /χ = n φ D /χ /s, where n φ D /χ and s are the actual number density of φ D and χ and entropy density of the universe are given by, , . (3. 2) The entropy density and Hubble parameter in terms of m φ D are where M r pl = 2.44 × 10 18 is the reduced Plank Mass. g * and g s * are the effective degree of freedom related to the energy and entropy density of the universe, respectively at temperature T = is the equilibrium number density of species i in comoving volume and is given by, with m i and g i are the mass and the internal degree of freedom for particle i, and K 2 the order-2 modified Bessel function of the second kind. The thermal average width Γ i of the species i is given by, The thermal average cross-section is given by [39], (3.6) where s is the centre of mass energy, and p ij (p kl ) are initial(final) centre of mass momentum. Finally, the relic density of the DM state χ is given by, where Ω T F I χ h 2 is the relic density obtained from the thermal freeze-in mechanism and Ω F O φ D h 2 is the abundance of φ D at the decoupling epoch T d . In the above, the second term represents the super-wimp contribution to the DM relic abundance, which occurs due to the late decay of φ D . The analytical expression for Ω T F I χ h 2 for the production of DM χ through the decay of φ D is given by, where g φ D is the internal degree of freedom for φ D . For the analysis, we consider the following mass spectra, M N = 50 GeV, M Z BL = 7 TeV, M H2 = 500 GeV, m φ D = 100 GeV (unless mentioned otherwise), and the gauge coupling g BL = 0.9. The choice of M Z BL and g BL is consistent with the constraint from CMS and ATLAS searches [40,41]. The right-hand side of Eq. 3.1 takes into account all possible number changing processes of φ D to B/SM states as well as its decay.
• This is to note, the depletion rate of φ D via Z BL mediated annihilation processes, i.e., φ † D φ D → Z BL → N N,f f, H 2 H 2 is suppressed due to large Z BL mass. Such processes decoupled from the cosmic soup much earlier compared to the annihilation of φ D through contact interactions and processes mediated via B/SM Higgs(H 1 , • The annihilation of φ D through contact interactions and s-channel mediated process via B/SM Higgs keeps the φ D in the thermal bath for a longer time. When the respective interaction rate becomes less than the universe's expansion rate, φ D decouples from the thermal bath.
• The depletion rate of φ D via processes that are dependent on the dark sector Yukawa coupling Y Dχ , such as χφ D → N Z BL , N φ D → χZ BL , etc are highly suppressed because of small coupling strength Y Dχ and negligible abundance of χ at an early epoch.
• The decay of φ D through φ D → χν process happens because of the active and sterile neutrino mixing which takes place after electroweak symmetry breaking (EWSB). The Heaviside step function, θ(x − x ew ) ensures that decrease in the number density of φ D via φ D → χν happens only after EWSB.
In Eq. 3.2, the right-hand side contains all relevant processes to study the evolution of χ. As discussed earlier, production of χ is dominated by the decay process, i.e., φ D → χN compared to the annihilation of the bath particles. Based on the primary production mechanism of χ, there are two different scenarios.
• Scenario-I : The DM is primarily produced via the thermal freeze-in mechanism. This corresponds to the case where φ D stays in the thermal bath for a more extended period of time owing to more considerable coupling strength with the bath particles. This tends to reduce its number density significantly. Thus its late decay gives negligible contribution to χ number density. Therefore, the correct relic density of χ is mostly obtained from the thermal freeze-in mechanism. This is illustrated in the left panel of Fig. 3. • Scenario-II : The DM χ is primarily produced from the decay of φ D after a freeze out, referred to as super-wimp mechanism. In this scenario, the number density of χ increases gradually through the thermal freeze-in mechanism, and at a later epoch, its number density increases significantly from the late decay of φ D . This is to note that this scenario is possible to realise if φ D decouples much earlier from the thermal bath due to suppressed interaction which leads to a large φ D abundance at the time of decoupling. φ D later decays completely to χ, significantly increasing the DM number density. We illustrate this schematically in the right panel of Fig. 3.
Before focusing on the main study of this work, we want to bring attention to the readers that φ D must decay before the big bang nucleosynthesis (BBN). The decay of φ D adds relativistic species to the thermal bath, which may alter the standard BBN scenario, and hence will be severely constrained. To avoid such conflict, we demand that lifetime of φ D must be less than 1 sec which puts a lower bound on dark sector Yukawa coupling Y Dχ , via which φ D decays. In Fig. 4, we show the lifetime contour of the φ D state for each of the given Yukawa coupling Y Dχ , where we vary the mass m χ of the DM and the mass of the parent particle m φ D . For each of the chosen Y Dχ values, the shaded region is allowed from BBN. We note that one of the products of the φ D decay is N , where N further decays to the SM particles and adds relativistic species to the thermal bath. The decay length of RHN depends on the active-sterile neutrino mixing, which depends on the light neutrino masses and PMNS mixing angles. We provide the expressions for the decay width of φ D and N in the appendix. For M N > m W ± , m Z , N decays dominantly through 2 body processes, such as, N → l ± W ∓ and N → Zν. For M N < m W ± , m Z which we consider in this study, RHN decays to the three SM fermions through off-shell W , and Z gauge bosons. For our choice of mass parameter M N , we have checked that the decay of N to SM states takes place instantaneously and is not constrained from BBN.  One can additionally set a upper bound on dark sector Yukawa coupling Y Dχ by demanding that φ D decays to χ after φ D freezes-out (at T ∼ m φ D /25). This requirement implies, In this work we consider such values of Y Dχ , so that both the BBN and the above mentioned constrained are satisfied. The entropy of the universe increases after the decay of φ D . The increase in entropy can be approximated [42] as, Due to the small value of Y Dχ , φ D decays in the late epoch of the universe when the universe temperature is around MeV. This leads to negligible amount of increase in entropy density in the universe, otherwise such increase in entropy can dilute the DM produced in an early epoch.

φ D abundance
Since a significantly large number of χ is produced from the late decay of φ D ; therefore the abundance of φ D at the time of decoupling plays a vital role in determining the correct DM relic density.
A large abundance of φ D will contribute significantly in the χ production via φ D → χN late decay.
As discussed before, φ D was in equilibrium with the SM particles because of the large scalar quartic couplings λ SD , λ Dh and SM-BSM Higgs mixing angle sin θ. The most dominant annihilation channels for φ D , φ † D φ D → W + W − , ZZ, N N mediated via SM-like Higgs state H 1 , for which we provide the analytic expressions of the cross-section in the appendix, and show the variation of thermal average cross-section σv F O at T d in Fig. 6(b). The H 2 mediated contribution is relatively smaller due to heavy propagator suppression, except the region of the H 2 resonance. In Fig. 5, we show the scatter plots in the λ SD and λ Dh plane for the three different values of mixing angle sin θ, for which .0 (for the bottom plot). The φ D abundance at T d is mostly governed by the couplings λ SD , λ Dh and sin θ. Comparing the horizontal and vertical bands between different panels, for a fixed value of sin θ, decreasing λ SD and λ Dh will lead to a higher Ω F O φ D h 2 . This typically occurs, as with the decrease in the relevant quartic coupling, the interaction rate of φ D decreases, resulting in an early freeze-out of φ D , which subsequently gives a larger φ D abundance. On the other hand, in each of these three panels, for a fixed value of λ Dh , as we decrease sin θ from 0.1 to 0.01 and further, a larger λ SD coupling is required to compensate the decrease in the interaction strength and to maintain Ω F O φ D h 2 in the given range. The cone-shaped region in each of the plots represents a cancellation in the φ † D φ D H 1 vertex that we will discuss later. Due to a relative suppression in the vertex factor, the interaction rate decreases, thereby leading to a higher value of in the given range, hence a larger value of couplings λ SD and λ Dh , and a large thermal averaged cross-section σv are required. To further explore the effect of an early and late decoupling of φ D on DM number density, we consider two benchmark points which are as follows, In Fig. 6(a), we show the variation of φ D abundance at T d with the change in φ D mass for these two above mentioned benchmark points. The red line in Fig. 6(a) corresponds to our first benchmark point, where φ D stays in thermal equilibrium for a longer period due to a large interaction rate with the bath particles which happen because of a large sin θ and λ SD . In Fig. 6 In each of these three panels, the green line corresponds to the second benchmark point, for which interaction of φ D is suppressed owing to a smaller coupling λ SD and mixing angle sin θ. The variation of Ω F O φ D h 2 and σv F O with the change in φ D mass is similar to the first benchmark point. It is important to note that due to suppressed interaction, φ D decouples from the thermal bath much earlier, leading to a large relic density of φ D . This is also reflected in Fig. 6(c), which shows the variation of the freeze-out temperature of φ D with its mass for these two benchmark points. As it is evident, the freeze-out temperature is relatively smaller for the first benchmark point, and hence freeze-out of φ D occurs at a later epoch. The sudden dip in the freeze-out temperature around φ D ∼ 60 GeV and 250 GeV occur because of s-channel resonance mediated via H 1 and H 2 . As the σv F O increases significantly in this region, this enables φ D to    remain in a thermal bath for a long time. The first benchmark point corresponds to Scenario-I, and the second benchmark point corresponds to Scenario-II, discussed earlier.

Vertex
Vertex Factor In Fig. 7(a), 7(b) and 7(c), we show the variation of φ D abundance at T d with the parameters λ SD , λ Dh and sin θ. We re-emphasize, the dominant annihilation mode for   contribute negligibly in φ D annihilation, as they are suppressed by the small mass of the final state fermions. In Fig. 7(a), we show the variation of Ω F O φ D h 2 with λ Dh while keeping λ SD and sin θ fixed However, for φ † D φD → N N , H1 mediated process is suppressed due to sin θ.
to few sets of values. As we can see from the green line which corresponds to λ SD = 10 −3 and sin θ = 10 −1 , Ω F O φ D h 2 remains independent of λ Dh in between 10 −5 to 10 −4 . This occurs as the effective vertex factor involving φ † D φ D H 1 for such a small value of λ Dh is governed by λ SD and sin θ rather than λ Dh . This can be understood from the expression of the vertex factors, which we provide in Table. 2. We also provide the vertex factors of different H 1,2 interactions with the SM fermions and gauge bosons in Table. 3. In between 10 −5 < λ Dh < 10 −4 , dominant annihilation modes are φ † D φ D → W + W − , ZZ mediated by SM Higgs boson H 1 . As λ Dh increases, effective vertex factor involving φ D φ † D H 1 decreases due to a relative cancellation between different terms in the respective expression (see Table. 2), leading to a suppressed annihilation rate. This results in an increase in the . This occurs as the process φ † D φ D → N N mediated by H 2 dominates as λ H2 becomes larger than λ H1 and also the process is not suppressed by scalar mixing angle. As λ Dh further increases, λ H1 increases which makes the annihilation rate of φ D large. This results in a decrease in Ω F O φ D h 2 as λ Dh in between 2×10 −3 < λ Dh < 10 −1 . The magenta and the yellow line follow the same characteristics as the green line. The difference in the abundance of φ D for different lines arises due to different choices of λ SD and sin θ. For a very large λ Dh , the green, yellow and magenta lines merge as λ H1 is governed by the λ Dh only. The red line corresponds to a much larger value of λ SD = 10 −1 and sin θ = 0.3. For this chosen parameter, λ H1 is governed by λ SD only. This makes the relic abundance of φ D at T d independent of λ Dh .
In Fig. 7(b), we show the variation of Ω F O φ D h 2 w.r.t the coupling λ SD . For the chosen parameters corresponding to the red, green and blue lines, Ω F O φ D h 2 is independent of λ SD for minimal coupling. In this scenario, λ H1 is governed by λ Dh only. As for both the green and red lines, the chosen value of λ Dh is the same and cos θ ≈ 1, therefore both of these lines give similar contributions to the Ω F O φ D h 2 for small λ SD . The blue line corresponds to a relatively smaller value of λ Dh compared to the green and red line, because of which Ω F O φ D h 2 is larger for blue line due to low interaction rate, in the region where λ SD is very small. As λ SD increases, the second term in the λ H1 becomes larger. This leads to cancellation in the respective vertex due to a difference in the sign between the first and second terms. Similar to Fig. 7(a), the suppression in the effective vertex leads to an increase in Ω F O φ D h 2 . In the cancellation region, the most dominant annihilation channel is φ † D φ D → N N . For the yellow line, λ Dh is considered to be negligible because of which λ H1 is governed by its second term only. Thus, as λ SD increases, the annihilation rate of φ D also increases which result in a decrease in Ω F O φ D h 2 . In Fig. 7(c), we show the variation of Ω F O φ D h 2 w.r.t sin θ. Similar to Fig. 7(a) and Fig. 7(b), for very small sin θ, few of the lines represent similar values of φ D abundance, as the interaction rate is totally governed by λ Dh cos θ combination in the λ H1 . The cancellation in λ H1 takes place only at a larger value of sin θ, for which Ω F O φ D h 2 increases significantly. For the chosen parameter corresponding to the red line, λ SD is much larger than λ Dh . Therefore, in this region, the annihilation rate is governed by the second term of λ H1 , because of which Ω F O φ D h 2 decreases for larger values of sin θ due to an increase in the annihilation rate. and Ω χ h 2 , we consider benchmark point 1. This is evident from Fig. 8(a) that φ D stays in the thermal bath for a significantly longer time owing to a large λ SD and sin θ. Due to this, the abundance of φ D is reduced significantly before it freezes out, and therefore, its contribution to the production of χ via late decay is negligible. For this figure, we consider a Y Dχ which is in agreement with the BBN constraint as well as Eq. 3.9. The lifetime of φ D increases with the increase in χ mass, thereby leading to the differences that can be seen from yellow, blue, green and red lines in the plot. In Fig. 8(b), we show the thermal freeze-in production of χ, production of χ from the out of equilibrium decay of φ D , and the relic abundance of χ including both the contributions. At a very early epoch, the abundance of χ was vanishingly small due to suppressed interaction of χ with bath particles owing to small coupling strength Y Dχ ∼ 10 −12 . For our chosen parameters, the production of χ is, however, primarily governed by the decay of φ D to χN state. The thermal freeze-in production of χ ceases as soon as the temperature of the thermal bath becomes less than the mass of φ D . In evaluating the relic abundance of χ, we have neglected the inverse decay process χN → φ D , as due to a very small abundance of χ, inverse decay is entirely negligible. Due to this, the freeze-in temperature of χ in our analysis is independent of abundance of χ but rather depends on m φ D . The abundance of χ can further be enhanced through the out of equilibrium decay of φ D . However, in this scenario, non-thermal production of χ from the late decay of φ D is tiny, as has been shown in Fig. 8. Therefore, the total production of χ, in this case, is determined by the thermal freeze-in mechanism. It is also important to highlight that the production of χ through the late decay of φ D increases as the mass of χ increases, as can be understood from Eq. 3.7. And it is also evident from the lower right side of Fig. 8(b), where we show the non-thermal contribution to relic abundance of Ω χ h 2 for different χ masses. Contrary to the previous scenario Scenario-I, Fig. 9 captures all the details about Scenario-II, where late decay of φ D contributes significantly in the production of χ. As it is evident from Fig. 9(a) that φ D decouples from the thermal bath much earlier due to suppressed interaction with the bath particles. The abundance of φ D at T d is therefore much larger than the φ D abundance for Scenario-I, and its out-of-equilibrium decay can contribute significantly to the χ abundance. This has been shown in Fig. 9(b), where we show the evolution of χ abundance. At high temperatures, the thermal freeze-in mechanism governs the production of χ. Similar to the previous scenario, at this very early epoch, the dominant freeze-in production mode of χ is from φ D decay. It is important to highlight that as the mass of χ increases from m χ = 50 GeV to m χ = 75 GeV, the thermal freeze-in contribution decreases instead of increase. It is due to the fact that phase space suppression for φ D → χN increases as mass of χ increases from m χ = 50 GeV to m χ = 75 GeV . At a later epoch, out of equilibrium decay of φ D starts to contribute to the production of χ. As one can see, the out of equilibrium decay of φ D alone can overproduce the χ. For our choice of parameters, the DM relic abundance is satisfied if m χ = 6.5 GeV. In Fig. 10, we show the effect of dark sector coupling Y Dχ on the production of φ D and χ. As one sees from Fig. 10(a) that change in Y Dχ only changes the lifetime of φ D . Owing to a small value of Y Dχ , it does not have any effect in the freeze-out processes of φ D . In Fig. 10(b), we show that the thermal freeze-in production of χ increases as Y Dχ increases. Since φ D abundance is independent of Y Dχ coupling, hence, production of χ from late decay of φ D is also independent of dark sector coupling Y Dχ . For the parameter choice, both the thermal freeze-in and non-thermal contributions are significant in the production of DM χ. Fig. 11 shows the variation of Ω F O φ D h 2 w.r.t x for different choices of m φ D . As one can see from the left panel ( Fig. 11(a)) that φ D yield increases as we vary m φ D = 100 GeV to m φ D = 150 GeV. However, for even larger m φ D values, such as, 200 and 250 GeV, the Ω F O φ D h 2 decreases as annihilation processes of φ D approaches s-channel resonances mediated via H 2 . This in turn effects the production of χ from the late decay of φ D . As one can see from Fig. 11(b), the thermal freeze-in production of χ at high temperature increases as the mass of φ D decreases. This occurs because T F I ∼ m φ D , and a lower T F I leads to higher production. At a later epoch, non-thermal production of χ from the late decay of φ D starts to contribute. The non-thermal production of χ increases as mass of φ D increases from 100 to 150 GeV, but later decreases with the increase in the mass of φ D .

3.4
Dependence of Y Dχ on model parameters In Fig. 12, we show the dependence of dark sector Yukawa coupling Y Dχ which governs the abundance of χ via thermal freeze-in production on parameters which determine the abundance of φ D at the time of decoupling. In Fig. 12(a), we show dependence of Y Dχ on m φ D for our two scenarios.  Figure 11. Left panel: Fig 11(a) shows the variation of Ω φ D h 2 w.r.t m φ D . Right panel: Fig 11(b) shows the variation of Ω χ h 2 w.r.t x.
• The thermal freeze-in dominated scenario, i.e., Scenario-I is represented by the red and green lines in the figure. In this scenario, φ D has a negligible abundance after freeze-out from the thermal bath, which is evident from the red line of Fig. 6(a). To show the dependence of Y Dχ on m φ D , we have considered two different masses of χ which are 10 and 20 GeV. The thermal freeze-in production ceases when the temperature of the thermal bath becomes less than m φ D . The freeze-in temperature drops as we consider a lighter φ D state. For the lower mass of φ D , production of χ, therefore, takes place for a longer time which in turn increase the abundance of χ significantly. To satisfy the correct relic density, we can decrease the production of χ by decreasing the dark sector Yukawa coupling Y Dχ . This behaviour is opposite for a large mass of φ D . The freeze-in production of χ suffers Boltzmann suppression at a much higher temperature compared to the case of the low mass of φ D . To compensate this effect, the production rate needs to be increased, which is done by increasing the magnitude of Y Dχ . It is important to highlight that as mass of φ D decreases from m φ D = 120 GeV to m φ D = 80 GeV, the magnitude of Y Dχ increases instead of decrease. It is due to fact that the phase space suppression for φ D → χN process increases as mass of φ D decreases from m φ D = 120 GeV to m φ D = 80 GeV. The behaviour of these two curves, even though similar, however the required value of Y Dχ is more prominent for the lower mass of χ compared to the higher mass of χ.
• The variation of the required Y Dχ for Scenario-II which satisfies the DM relic abundance is shown by the pink and blue lines in Fig.12(a). Ω F O φ D h 2 for Scenario-II is shown via the green coloured line in Fig. 6(a), which indicates Ω F O φ D h 2 ≥ 0.12 except the region near schannel resonance around m φ D = M H1 /2 ∼ 62.5 GeV and m φ D = M H2 /2 ∼ 250 GeV. The large abundance of φ D enhances the relic density of χ. The late decay of φ D producing χ is independent of the Yukawa coupling Y Dχ ; however, the thermal contribution depends on Y Dχ . Hence, depending upon the abundance of φ D , the dark sector Yukawa coupling Y Dχ needs to be tuned accordingly to satisfy the relic density constraint. It is important to highlight that the out-of-equilibrium decay of φ D in the production of χ can be so significant that to satisfy the correct relic abundance of χ, the thermal freeze-in production of χ is required to be small, which is possible to achieve for a small Y Dχ . However, note that the coupling Y Dχ can not be made arbitrarily small, as the BBN imposes a strong lower bound on Y Dχ . The pink and blue lines in Fig. 12(a) along which the DM relic abundance Ω χ h 2 = 0.12 show the variation of the required Y Dχ w.r.t φ D mass for this scenario. The pink line clearly shows that a smaller value of Y Dχ is required in order to satisfy the correct relic density for Scenario-II when compared to the red line, which corresponds to the thermal freeze-in dominated scenario of Scenario-I. Also, note that the pink and red lines merge for the mass of φ D greater than 200 GeV. For m φ D ∼ 250 GeV, due to the s-channel resonance, the abundance of φ D decreases significantly (see Fig. 6(a)). Therefore, a large thermal freeze-in contribution is required to satisfy the correct relic abundance, which in turn demands a larger value of the Yukawa coupling Y Dχ . The blue line traces the pink line in part of the parameter space, with the notable difference that the out of equilibrium decay of φ D producing χ is more dominant in between 110-190 GeV. The relic abundance of χ in this region becomes larger than the observed relic density due to the large contribution from the late decay of φ D . Hence, for no value of Y Dχ , the DM relic density constraint Ω χ h 2 = 0.12 is satisfied.
In Fig. 12(b), we show dependence of Y Dχ on λ Dh . The observations are listed as follows: • The red line in this figure represents Scenario-I, i.e., the thermal freeze-in dominated scenario.
In Fig. 7(a), the red line shows the variation of Ω F O φ D h 2 with λ Dh for λ SD = 10 −1 and sin θ = 0.3. Ω F O φ D h 2 is significantly small for all values of λ Dh . Therefore, out of equilibrium decay of φ D can not produce significant number of χ. Due to this, the correct relic density of χ is obtained only through thermal freeze-in production which depends on Y Dχ , and not on the coupling λ Dh . Therefore, the required coupling Y Dχ is independent of λ Dh .
• The variation of Y Dχ w.r.t the variation of λ Dh in Fig. 12(b) can be understood from Fig. 7(a).
For fixed value of m φ D and m χ , production of χ through out of equilibrium decay of φ D is proportional to Ω F O φ D h 2 . As we can see from the green and yellow lines in Fig. 7(a) that Ω F O φ D h 2 > 0.12 but remain constant in the range 10 −5 < λ Dh < 10 −4 . For λ Dh > 10 −4 , Ω F O φ D h 2 increases significantly as λ Dh increases. This sudden jump occurs due to the cancellation in λ H1 , described earlier. Ω F O φ D h 2 again falls much below 0.12 as λ Dh increases further. For a large φ D abundance, the out-of-equilibrium contribution from φ D → χN will be substantial. In Fig. 12(b), for λ Dh > 3 × 10 −3 , due to very suppressed φ D abundance, out of equilibrium decay contribution is small, and the thermal freeze-in contribution alone satisfies the DM abundance. This is represented via the red, yellow and green lines which merge near λ Dh ∼ 3 × 10 −3 . For small value of λ Dh < 10 −4 , both the out of equilibrium decay of φ D as well as thermal freeze-in production contribute substantially to the relic abundance of χ. In this region, due to the presence of a finite out-of-equilibrium decay contribution, the required value of Y Dχ to satisfy correct DM relic density is typically less when compared to the only thermal freeze-in dominated scenario, represented via the red line. In between 5×10 −4 < λ Dh < 2×10 −3 , the φ D abundance is very large due to cancellation in λ H1 leading to a large out-of-equilibrium decay contribution which results in Ω χ h 2 > 0.12. Hence, this region is disallowed.
• The magenta line in this figure corresponds to λ SD = 10 −3 and sin θ = 10 −2 . The behaviour can again be understood by referring to the magenta line in Fig. 7(a). Ω F O φ D h 2 is much larger than 1.0 for λ Dh < 2 × 10 −3 . This leads to the overproduction of χ via out-of-equilibrium decay of φ D . Hence, the correct relic density of χ in this case is only obtained for λ Dh larger than 2 × 10 −3 .
In Fig. 12(c), we show dependency of Y Dχ on λ SD . The observations are listed as follows: • For the yellow curve in Fig. 12(c), the correct relic density of χ is possible to obtain for λ SD > 3 × 10 −4 . For λ SD < 3 × 10 −4 , abundance of φ D is significantly large (see the yellow curve of Fig. 7(b)), leading to the overproduction of χ via out of equilibrium decay of φ D .
• We can see in Fig. 7(b) that red, green and blue lines have similar features. The cancellation in λ H1 takes place for different values of λ SD . As we can see the funnel shaped region in Fig.12(c) that Y Dχ decreases significantly in the region where cancellation in λ H1 is effective.
In Fig. 12(d), we show dependence of Y Dχ on sin θ. The observations are listed as follows: • The nature of the red, green and blue lines can be understood from Fig. 7(c). The required value of Y Dχ to satisfy correct DM relic density decreases significantly when φ D abundance gets enhanced. For each of these three lines, in the funnel shaped region, the φ D abundance becomes very large due to the cancellation in λ H1 . This larger φ D abundance leads to Ωh 2 χ > 0.12, which is ruled out. • In Fig.12(d), the red line represents a scenario, where both the thermal freeze-in production and out-of-equilibrium production of χ can contribute. For sin θ > 0.1, due to a smaller φ D abundance, mostly thermal freeze-in contribution dominate. Hence a larger value of Y Dχ is required to satisfy the DM relic abundance. In the blue, pink, and green lines, the effect of cancellation in the φ † D φ D H 1 vertex is clearly visible. For each of these lines, in the funnel shaped region, φ D abundance is very large, leading to an overproduction of χ. For small sin θ, both the thermal freeze-in and out-of-equilibrium decay can contribute significantly (see Fig. 7(c)).

Collider Prospects
This section focuses on the search for the BSM Higgs H 2 via its invisible decay, i.e., H 2 → φ † D φ D at the LHC. The partial decay width for this decay mode is determined by the couplings λ SD . The other production mode of φ D from Z BL is suppressed due to a very heavy Z BL . The possible decay mode of φ D is φ D → χN which is controlled by the coupling Y Dχ . As mentioned before we assume Y Dχ = O(10 −12 ) to realize the freeze-in production of the DM χ. As a result of this tiny coupling, φ D escapes the detector without leaving any visible footprint. However, its production can be confirmed by the observed imbalance in the transverse momentum. For collider analysis, we consider the mass of φ D to be m φ D = 100 GeV. Other parameters are set to M Z BL = 7 TeV, Parameter κ Z κ W κ t κ b κ ta κ g κ γ sin θ at 95%CL 0.46 0.45 0.65 0.69 0.65 0.61 0.42 Table 4. Upper limit on sin θ obtained from Higgs boson coupling modifiers, κ Z/W/t/b/τ /g/γ at 95% CL [43].
g BL = 0.9, and M N = 50 GeV, which we also consider for the DM study. We first discuss existing constraints on the model parameters from the collider experiment and then project the future sensitivity to probe the coupling λ SD at the HL-LHC.

LHC Constraints
We first discuss the different constraints applicable on the SM-BSM Higgs mixing angle θ, the mass of the BSM Higgs and the quartic coupling λ SD .

Measurement of Higgs signal strength and coupling constant modifiers:
The signal strength of SM Higgs decaying into two SM states a, b, such as W W * , ZZ * , τ τ, bb, µ + µ − is, The global signal strength of H 1 using 139 fb −1 data at LHC is measured as µ = 1.06±0.07 [43]. Due to the presence of the BSM Higgs, which mix with the SM like Higgs states, the standard couplings of the SM Higgs with W + W − , ZZ, τ τ and others will be modified. We adopt the constant coupling modifiers -κ framework, where κ's are defined as where λ xxh is the couplings of SM-like Higgs field H 1 with two SM fields in the model considered, and λ SM xxh is the respective coupling in the SM. We consider the ATLAS search [43], and translate the measurements of each measured κ's to the upper limit on the SM and BSM Higgs mixing angle sin θ. The results are shown in Table. 4. In our collider analysis for HL-LHC, we adopt a conservative approach and consider relatively smaller values of sin θ = 0.3, which agrees with the LHC constraints. SM Higgs decaying to long-lived particle (LLP): The theory under consideration predict several exotic decays of the SM Higgs boson, such as and its possible decay modes are N → ljj/νjj/llν/3ν. Therefore, N is a LLP undergoing displaced decays. The recent CMS search for displaced heavy neutral leptons limits the active-sterile mixing in the mass range 1−18 GeV [44] with the most tight constraint appears |V | 2 < 10 −7 for M N ∼ O(10) GeV. Our choice of RHN mass M N = 50 GeV and active sterile mixing V ∼ 10 −7 is beyond the range covered in this paper. There are other CMS and ATLAS searches for exotic decays of SM Higgs into two LLP states, which are instead applicable. The CMS and ATLAS have recently searched for exotic decays of the Higgs boson into LLP in the tracking system [45,46]. These searches are mainly sensitive to LLP with cτ = O(1 mm − 300 mm). Other displaced vertex searches in the tracking system that are also sensitive to Higgs decays to  Fig. 13(b): Constraints in sin θ − λ SD plane derived from the ATLAS search for heavy scalar resonance decaying to two Z bosons, pp → H 2 → ZZ → 4l [53].
LLP [47,48]. Our benchmark point is unconstrained from these searches owing to a very long lifetime of the RHN. The latest search for neutral LLP decaying into displaced jets in the ATLAS muon spectrometer [49][50][51] and in the CMS endcap muon detectors [52] are relevant for LLP with cτ ≥ O(1 m). The RHN mostly decays in the muon spectrometer for our benchmark mass point. This is to note that our model prediction of BR(H 1 → N N ) = 0.5% for M N = 50 GeV is consistent with the observed bound on BR of Higgs to LLP decay. Note that this constraint is given for the Higgs decaying to scalar LLP and for the two-body decay of the LLP. In reinterpreting this analysis for our scenario, we assume similar signal selection efficiency as given in [52].
Heavy Higgs searches: Other LHC searches [53][54][55][56][57] aimed at probing BSM Higgs via direct measurements can constrain our model. These are the searches to detect a heavy scalar resonance (H 2 ) decaying into various final states, such as W + W − /ZZ//H 1 H 1 . Among them the strongest limits comes from the multi-lepton search in the channel pp → H 2 → ZZ [53]. In the model under consideration, H 2 has an additional decay mode H 2 → φ † D φ D , which is governed by the coupling λ SD . For large value of the coupling λ SD this decay mode can be dominant over H 2 → H 1 H 1 /W + W − /ZZ. Fig. 13(a) shows the contours of BR(H 2 → φ † D φ D /W + W − /ZZ//H 1 H 1 /tt) in the M H2 -λ SD plane. The expressions for the respective partial decay widths are given in the appendix. For this plot we fix the scalar mixing angle, sin θ = 0.3. The gray shaded region represents grows with λ SD , BR(H 2 → ZZ) decreases and hence, the bound on scalar mixing angle sin θ becomes weaker for a fixed mass of H 2 . This is shown in Fig. 13(b) for three illustrative mass points of H 2 , M H2 = 350, 500, 1000 GeV. Here we translate the observed limit on σ(pp → H 2 → ZZ) from ATLAS search [53] into sin θ − λ SD plane. The shaded regions are disallowed for the respective values of M H2 . For a smaller value of λ SD < 10 −2 , for which H 2 → ZZ branching is significantly larger, a very tight constraint sin θ < 0.2 appears for M H2 < 500 GeV.
The CMS and ATLAS collaborations have also performed searches for Higgs boson decaying invisibly. For our parameter choice H 1 → φ † D φ D is closed. Recently ATLAS has searched for such invisible decay of Higgs through vector boson fusion (VBF) production channel and interpreted the result for a heavy scalar particle [58]. In Fig. 14, the black curve shows the observed bounds on the cross-section times branching ratio to invisible final states of the heavy Higgs from this ATLAS search. The blue-dashed curve represents the theory prediction for BR(H 2 → inv) = 1 and for the scalar mixing angle sin θ = 0.3. In deriving this, we consider a simplistic parton-level analysis with MadGraph5 and do not consider any specific cut-efficiencies. Our theory cross-section agrees with the observed limit in the entire mass range, which is evident from this plot. In the upcoming section, we examined the reach of HL-LHC to search for H 2 decaying invisibly through VBF production mode. The Feynman diagram for this process is shown in Fig. 15.

Search for
The VBF process is one of the most promising channels to search for the invisible decay of Higgs boson [59]. Recently the ATLAS [58] and CMS [60] collaborations have studied the SM Higgs decay to invisible particles and constrained such production processes. Invisible decay of the SM Higgs boson through the VBF channel has been studied for the Higgs portal models [61,62], for Inert-doublet model [63]. Below, we investigate the production of the BSM scalar H 2 via the VBF process and its subsequent decay to the invisible state φ D . Note that, owing to a very tiny coupling Y Dχ ∼ O(10 −12 ), φ D state decays outside the detector. After gluon-fusion, VBF is the dominant channel for Higgs bosons production at the LHC, characterised by the two highly energetic forward jets [64]. The two VBF jets are widely separated in pseudo-rapidity, lying in the opposite hemisphere of the detector. For the invisible decay H 2 → φ † D φ D , the signal is marked by a large transverse momentum imbalance. All these features allow us to discriminate between the signal and background. The dominant SM processes that mimic the signal are pp → Z(→ νν)jj and pp → W ± (→ ν ± )jj. The latter process contributes when the charged lepton is not detected. QCD multi-jet events with large missing transverse momentum (MET), arising from the mismeasurement of jet energy, can also imitate the VBF signal. A suppressed central jet activity accompanies the VBF signal. On the contrary, QCD jets are more central in the detector. Therefore, the central jet veto and a strong cut on MET could reduce QCD multi-jet events. Another potential background tt can be suppressed by vetoing b jets and leptons. showering and hadronization. We implement a cut-count analysis code in CheckMate [68,69], to calculate the signal and background cut efficiencies. CheckMate makes use of Delphes [70] for the simulation of detector effect, and Fastjet [71,72] for jet clustering. We use anti-k t jet clustering algorithm [73] with radius parameter, R = 0.4. We estimate the sensitivity of the invisible signature of the H 2 produced via VBF process at the pp-collider, pp → H 2 jj → φ † D φ D jj. Here φ D being a stable particle at the detector length scale gives rise to MET. Thus, the process under consideration leads to 2j + M ET signature at collider. Among the other SM processes that can fake the signal we simulate the two most dominant processes which are pp → Zjj → ννjj and pp → W ± jj → ν ± jj. We consider the HL-LHC for this study, which is planned to operate with √ s = 14 TeV and L = 3000/fb.
Although the signal consists of two jets at the parton level, additional jets can arise due to initial and final state radiation after the parton shower. Thus, we consider up to two extra jets in the final state to simulate backgrounds. During the generation of background events we demand transverse momentum (p T ) of the leading partons: p T (j 1,2 ) > 50 GeV, the pseudo-rapidity: |η j | < 5.0. We simulate merged sample with 2-4 jets for W + jets and Z + jets using the MLM matching scheme [74]. We consider the parameter xqcut=55 GeV which is the minimum jet measure (p T /k T ) between partons. Partonic cross-sections for backgrounds W + jets and Z + jets are 528.755 × 10 3 fb and 223.603 × 10 3 fb, respectively. Next to leading order (NLO) QCD corrections for W + jets and Z + jets are given in [75], which are negative, and the corresponding K-factors are 0.87 and 0.905 , respectively. We perform the simulation with a harder p T cut on jets compared to that in the mentioned reference. As, K-factor depends on the kinematic cuts we do not normalise the cross-sections to NLO. For signal we also consider LO cross-section. For M H2 = (350 − 1000) GeV, K-factor varies in the range 1.018 − 0.979 [76]. The sensitivity can be improved for reduced background cross-sections and enhanced signal cross-section. Fig. 16 shows the normalised distributions of some of the kinematic variables for both signal and background which motivate to design the selection cuts. Fig. 16(a) corresponds to the distributions for transverse momentum of the leading jet (p T (j 1 )) which shows that the background events possess comparatively harder jets than signal events due to the generation level p T cut. Fig. 16(b) and Fig. 16(c) represent the difference in pseudo-rapidity (|∆η(j 1 , j 2 )|) and invariant mass (M (j 1 j 2 )) of the leading and sub-leading jets. Majority of the signal events hold higher |∆η(j 1 , j 2 )| with respect to the background events and hence, larger M (j 1 j 2 ) as both variables are connected with the equation M (j 1 j 2 ) p T (j 1 ) p T (j 2 ) e ∆η(j1,j2) 4 .
Event selection and Results:-We closely follow the cuts used in ATLAS search [58], which are as follows  Figure 16. Normalised Distribution of transverse momentum of the leading jet ( Fig. 16(a)), difference in pseudo-rapidity ( Fig. 16(b)) and invariant mass of the leading and sub-leading jets ( Fig. 16(c)), for both signal and background events. The distributions are without any selection cuts. .
• We veto events with more than one b-tagged jets.
• We demand missing transverse momentum E miss T ≥ 150 GeV.
In Table. 5, we present the cut efficiencies of the cuts mentioned above for signal and the SM background. We find the lepton veto to be most effective in reducing the W ± + jets events. Finally, after the cuts on ∆η(j 1 , j 2 ) and M (j 1 j 2 ) a significant fraction of both W ± + jets and Z + jets events are cut down. The effect of these two cuts is similar as both variables are related. In the last row, we write the required signal cross-section (in fb) before applying selection cuts to obtain 2σ significance for L = 3000/fb luminosity. This is calculated from the relation is the initial signal (background) cross-section, s and b is the corresponding cut efficiency, and n σ is the significance.
In Fig. 17, we show the HL-LHC prediction for invisible signature of H 2 in M H2 − λ SD plane assuming sin θ = 0.3. The signal significance is calculated for 3000/fb luminosity. The brown solid (brown dashed-dot) line indicates to 2σ (5σ) sensitivity. As per our results, 2σ significance can be achieved upto M H2 800 GeV for λ SD in between 10 −2 and 1. We also highlighted the region excluded from the ATLAS search for H 2 → ZZ [53] for 2σ significance with L = 3000/fb 13.831 fb 11.905 fb 10.802 fb -- Table 5. Cumulative cut efficiencies for Signal: pp → H 2 jj → 2j + M ET and the SM background: W + jets and Z + jets. The numbers in last row are the required signal cross-sections to obtain 2σ significance using 3000/fb luminosity.
area enclosed by it corresponds to Γ/M H2 < 0.1. Note that to obtain the results for the invisible search we take into account the width effect. The sensitivity of the visible signature of H 2 has been analysed in ref. [77]. Here the authors have considered the production of H 2 via gluon fusion process and subsequent decay to di-Higgs, pp → H 2 → H 1 H 1 . They have studied various final states depending on the decay of SM Higgs H 1 and have shown that 2b + 2γ and 4b signatures are more sensitive compared to others. In Fig. 17, we shows the HL-LHC predictions for di-Higgs channel in 2b + 2γ and 4b final state obtained from the ref. [77]. Here we assume sin θ = 0.3. For 2b + 2γ signature, regions enclosed by red solid and dashed-dot curves represent 2σ and 5σ sensitivity, respectively. Similarly blue curves show 2σ and 3σ sensitivity for 4b final state.

Conclusion
We analyse the thermal freeze-in and non-thermal freeze-in production of DM in an extended, gauged B − L model where dark sector fermion χ serves as the DM candidate. In this work, we have a secluded dark sector containing feeble interacting DM candidate χ and a complex scalar field φ D charged under B − L symmetry. The DM fails to thermalise with the surrounding plasma due to suppressed interaction with other bath particles owing to small coupling Y Dχ . At an early epoch, it is primarily produced via the decay of φ D via thermal freeze-in mechanism and through the late decay of φ D via the non-thermal freeze-in mechanism. Contrary to the DM field χ, the dark sector scalar field φ D thermalises with the bath particles due to large gauge coupling as well as sizeable interactions with SM and BSM Higgs states. The correct abundance of φ D at decoupling is therefore obtained via the freeze-out mechanism. The annihilation of φ D mediated via Z BL , i.e., φ † D φ D → Z BL → ff etc decouples from the thermal bath at an earlier epoch compared to the processes mediated via SM and BSM Higgs. The abundance of φ D is hence primarily governed by the interplay of scalar quartic couplings λ SD , λ Dh and SM-BSM Higgs mixing angle sin θ. For our choice of φ D and χ masses, φ D decaying to χ and RHN is kinematically open, and this is the primary production process for DM. We subdivide the discussion into two scenarios, Scenario-I, where the thermal freeze-in via φ D → χN contribution is significant, and Scenario-II, where the late decay of φ D → χN is primarily responsible in satisfying the DM relic abundance. This is to note that the late production of DM through out-of-equilibrium decay of φ D will depend on the abundance of φ D at decoupling obtained through the freeze-out mechanism. Therefore, for suppressed interaction of φ D with the bath particles leading to high abundance of φ D at decoupling significantly enhances the production of χ. We also study mixed scenarios, where the correct relic abundance of χ is governed via both the thermal and non-thermal freeze-in mechanism. We majorly focus our discussion around two benchmark points, 1 and 2. We find that • Due to relatively large SM-BSM Higgs mixing angle, as well as due to the choice of large quartic couplings λ SD , λ Dh , φ D decouples much later in Scenario I, and hence its late decay contributes negligibly to the DM abundance. Here, the primary production of χ occurs due to thermal freeze-in.
• For benchmark 2 (Scenario II), the small SM-BSM Higgs mixing and choice of λ SD and λ Dh enables φ D to be out-of-equilibrium in an earlier epoch leading to a larger φ D abundance. In this case, the χ abundance can be built-up primarily from the late decay of φ D .
• We find that for Scenario-I, which corresponds to benchmark 1, the dark sector scalar φ D can be produced at a pp collider. This occurs because the SM-BSM Higgs mixing angle is relatively larger to realise this scenario. The φ D is produced from BSM Higgs decay in the VBF channel.
In addition to the detailed DM analysis, we also evaluate the prospect of detection of this model at the HL-LHC. In particular, in Section. 4, we investigate the possibility of probing the coupling λ SD of φ D with the heavy scalar H 2 at the HL-LHC. This coupling enables an extra decay mode of H 2 to a pair of φ D . When this decay mode becomes dominant, the existing bounds on the mass of H 2 and the scalar mixing angle weaken. To study H 2 → φ † D φ D , we consider the production of H 2 from the VBF process, characterised by two forward jets with a large pseudo-rapidity gap. In our case, φ D is stable over the detector length scale resulting in an extra distinguishing feature, the missing transverse momentum. For a fix mass of φ D , we present the 5σ discovery and 2σ exclusion contours in the M H2 − λ SD plane. Following a simple cut-count analysis we show that 5σ sensitivity can be obtained for M H2 (350 − 500) GeV and λ SD in between 0.03 and 0.35. Similarly, 2σ exclusion limit can be placed for the mass range (350 − 800) GeV and for λ SD (10 −2 − 1). where in the above,λ(x 2 , y 2 , z 2 ) = x 2 + y 2 + z 2 − 2xy − 2yz − 2zx is a Kälen function.

A.2 Decay width of N
The two body decay width of N R when m N is larger than m W , m Z and m H1 are as follows, • Γ(N → χφ) = For M N < m W ± , m Z , it decays to the three SM fermions through off-shell W , and Z gauge bosons. The three body decay widths are as follows, Here, I(x u , x d , x l ) = 12 Here x = . The values of C f 1 and C f 2 are given in [78] • In the above, n c is the color charge and is 1 for leptons and 3 for quarks.
Similarly, the expression of cross-section for the relevant processes in χ production are as follows,