Singlet-Doublet Fermionic Dark Matter and Gravitational Wave in Two Higgs Doublet Extension of the Standard Model

We present a study of singlet-doublet vector-like leptonic dark matter (DM) in the framework of two Higgs doublet model (2HDM), where the dark sector is comprised of one doublet and one singlet vectorlike fermions (VLFs). The DM, that arises as an admixture of the neutral components of the VLFs, is stabilized by an imposed discrete symmetry $\mathcal{Z}_2^{'}$ . We test the viability of the model in the light of observations from PLANCK and recent limits on spin-independent direct detection experiments, and search for its possible collider signals. In addition, we also look for the stochastic gravitational wave (GW) signatures resulting from strong first order phase transition due to the presence of the second Higgs doublet. The model thus offers a viable parameter space for a stable DM candidate that can be probed from direct search, collider and GW experiments.


I. INTRODUCTION
Despite strong evidence of the existence of the dark matter (DM) from several astrophysical and cosmological observations like rotation curves of spiral galaxies [1,2], the bullet cluster [3], gravitational lensing [4] etc., particle nature of dark matter (DM) is unknown till date. While Planck [5] results claim that nearly 26.5% of matter in the Universe is indeed DM, the Standard Model (SM) of particle physics is incapable of accounting for a viable DM candidate. Interestingly, if DM interactions with the SM particles are similar to those of electroweak interactions, and the particle DM has a mass around the electroweak scale, then such a DM can be thermally produced in the early universe, followed by its freeze-out, leaving a thermal relic very close to the observed DM abundance (Ωh 2 ∼ 0.12). This remarkable coincidence is often referred to as the weakly interacting massive particle (WIMP) miracle [6]. Although WIMP remains as an elusive DM candidate (for a recent review on the status of WIMP see [7]) as it indicates new physics signature around TeV-scale, but the non-observation of any excess both at the colliders and at the DM scattering experiments such as LUX [8], PandaX-II [9,10] and XENON1T [11,12] etc. compels us to strive for either an alternative to WIMP-paradigm [13][14][15][16] or to come up with some search strategies for DM detection other than the usual scattering experiments (for an overview see [17]).
Motivated from these, in this work we propose a particle DM model by extending the SM with a second Higgs doublet, and adding one vector-like lepton doublet and one vector-like lepton singlet. We consider a Type-I two Higgs doublet model (2HDM) where one of the Higgs doublets is odd under a discrete Z 2 symmetry. Such an imposition of a discrete symmetry is very commonplace in context with 2HDM as this necessarily prevents the appearance of tree-level flavour changing neutral current (FCNC) via Higgs [18][19][20][21][22][23]. The newly added fermions are assumed to be odd under a second Z 2 symmetry. All the SM particles are even under both of the discrete symmetries Z 2 and Z 2 . The imposition of two different discrete symmetries ensures a first order phase transition (FOPT) without hampering the stability of the DM. Under these circumstances the DM emerges as the lightest particle odd under Z 2 due to an admixture of the neutral component of the doublet and the singlet. Studies of minimal singlet-doublet vectorlike DM in SM has been exhaustively performed in the literature [24][25][26][27][28][29][30][31][32][33][34][35][36][37][38][39][40][41][42]. As it is understandable, a purely singlet vector-like fermion (VLF) DM does not have any renormalizable portal interaction with the SM to obtain the observed thermal relic abundance. The purely doublet VLF, on the other hand, annihilates too much to the SM due to its electroweak gauge interactions, thus making it under abundant unless the mass is > ∼ TeV. A purely doublet VLF also faces stringent constraints from DM direct detection experiments because of large scattering cross-section mediated by Z boson. In the present model, as we shall see, due to the presence of the second Higgs doublet the bound from direct search is less stringent. This is possible due to some destructive interference 1 between the scalar mediated direct search diagrams. Even in the absence of such destructive interference, a small direct search cross-section is still conceivable due to the suppression coming from the heavy Higgs mass and small scalar mixing. This provides some freedom of choosing a moderate sin θ. However, as the Z-mediated DM-nucleon scattering is still present, hence the constraint is not completely alleviated, and that confines the singlet-doublet mixing to some extent. This, in turn, affects the collider signature for this model. Here we would like to mention that our model is different from the one present in [40] where the DM particles do not couple to SM Z boson due to its Majorana nature.
The origin of the baryon asymmetry of the Universe (BAU) is another long-standing puzzle of particle physics. The electroweak baryogenesis (EWBG) is a possible way to account for the BAU exploiting the three Sakharov conditions [47]. However, it is not possible to have a successful EWBG within the SM paradigm as the SM neither provides sufficient CP-violation or strong firstorder phase transition (SFOPT) [48][49][50][51]. Therefore, a successful EWBG invokes new physics at the electroweak scale that can be obtained via an extended scalar sector. The two Higgs doublet model (2HDM) is a very well motivated non-supersymmetric extension of SM where the scalar sector of the SM is augmented with additional Higgs doublet giving rise to a plethora of different phenomenological implications [18][19][20][21][22][23][52][53][54][55][56]. In context with electroweak phase transition (EWPT), 2HDM has been extensively studied both in the CP-conserving case [57][58][59] and CP-violating scenario [60][61][62][63]. It was also shown that 2HDM framework is capable of generating a SFOPT [64,65].
The GW signature as a complimentary search strategy in context with DM models has already been studied both in case of freeze-out and freeze-in [67,[87][88][89][90][91][92][93][94](for a review on GW probes of DM see [95]). We have shown, within a consistent framework, our model is also capable of providing a detectable GW signal by satisfying all stringent DM, collider and other theoretical constraints.
The paper is organised as follows: in Sec. II we have introduced the particle content of the model 1 The presence of "bind spot" in direct search cross-section for a 2HDM model has been studied in [40,[43][44][45][46] along with the necessary interaction terms, in Sec. III we have discussed the constraints on the   model parameters arising due to tree-level unitarity, precision observables and collider bounds,   we next move on to Sec. IV where we illustrate the parameter space satisfying relic abundance and direct detection bounds from which we choose a couple of benchmark points to perform the collider analysis in Sec. V, then in Sec. VI we detail the generation of gravitational wave due to SFOPT Table I. In this set-up the Lagrangian for the model can be written as: where L f is the Lagrangian for the VLFs, L s involves the SM doublet and the additional Higgs doublet, and L yuk contains the Yukawa interaction terms.
The interaction Lagrangian for the VLFs reads: where D µ is the covariant derivative under SU (2) × U (1): where g and g are the gauge couplings corresponding to SU (2) and U (1) Y and a = 1, 2, 3 are the indices for the generators of SU (2). W µ and B µ are the gauge bosons corresponding to SM SU (2) and U (1) Y gauge groups respectively.
Lagrangian of the scalar sector involving SM Higgs doublet (H) and the new Higgs doublet (Φ 2 ) can be written as: The model with such a modified scalar sector thus resembles with the standard two Higgs doublet model (2HDM) of Type-I [18,[20][21][22], where all SM fermions have Yukawa interactions with only one of the doublets e.g., Φ 2 . The most general renormalizable scalar potential can then be written as: As we shall see, the coefficient m 12 of the softly Z 2 -breaking term plays the pivotal role in deciding the nature of the phase transition. Finally, the charge assignment allows to write a Yukawa interaction [32,37,41]: where Y is the Yukawa coupling between the VLFs and SM Higgs and Φ 2 = iσ 2 Φ * 2 .

B. Mixing in the scalar sector
We parametrize the scalar doublets as: for i = 1, 2. After spontaneous symmetry breaking (SSB), doublet Higgs fields acquires vacuum GeV. The ratio of VEVs is given as tan β = v 2 v 1 . The physical states are obtained by diagonalizing the charged and neutral scalar mass matrices. There are then altogether eight mass eigenstates, three of which become the longitudinal components of the W ± and Z gauge bosons. Of the remaining five, there is one charged scalars H ± , two neutral CP even scalars h, H and one neutral pseudoscalar A. The mixing between CP even scalars is denoted by an angle α . For our analysis we shall follow the where we denote s α = sin α, c α = cos α and similarly s β = sin β and c β = cos β.

C. Mixing in the VLF sector
After EWSB, the neutral components of the doublet (ψ 0 ) and singlet (χ 0 ) mix via the Yukawa interaction (Eq. 6). The mass matrix can be diagonalized in the usual way using a 2 × 2 orthogonal rotation matrix to find the masses in the physical basis (ψ 1 , ψ 2 ) T : where the non-diagonal terms are present due to Eq. 6, and the rotation matrix is given by The mixing angle is related to the masses in the weak (flavour) basis: The physical eigenstates (in mass basis) are, therefore, a linear superposition of the neutral weak eigenstates. These can be expressed in terms of the mixing angles as: The lightest electromagnetic charge neutral Z 2 odd particle is a viable DM candidate of this model. From now on we shall refer ψ 1 as the lightest stable particle (LSP) of the model. In the small mixing limit, the charged component of the VLF doublet ψ ± acquires a mass as: From Eq. 14, we see that the VLF Yukawa is related to the mass difference between two physical eigenstates and is no more an independent parameter: Therefore one can have three new parameters: {m ψ 1 , ∆m, sin θ} from DM phenomenology apart from 2HDM parameters. These three parameters will play the key role in determining the relic abundance of the DM, also deciding the fate of the model in direct and collider searches.

III. CONSTRAINTS ON THE MODEL PARAMETERS
In this section we would like to summarize constraints on the masses, mixings and couplings arising in the model due to theoretical and experimental bounds. We are particularly interested in the choice of tan β and heavy scalar masses in the 2HDM sector in order to have a SFOPT, while the free parameters appearing in the VLF sector is mostly constrained by oblique parameters and later from DM phenomenology.

Vacuum Stability
Stability of the 2HDM potential is ensured by the following conditions [21,22,97], These conditions have been shown to be necessary and sufficient [98] to ensure that the scalar potential is bounded from below [99].

Perturbativity
Tree-level unitarity imposes bounds on the size of the quartic couplings λ i or various combinations of them. The quartic couplings and the Yukawa couplings appearing in the theory need to satisfy [21,22,97]: in order to remain within the perturbative limit. Here λ i = λ, λ 1,2,3,4,5 . Here we would like to mention that absolute stability of the vacuum at tree-level puts a bound on tan β depending on the choice of m H as derived in [100]. For m H 300 − 400 GeV this typically allows 1 < ∼ tan β < ∼ 30.

Electroweak precision observables (EWPO)
The splitting between the heavy scalar masses is constrained by the oblique electroweak Tparameter whose expression in the alignment limit is given by [21,97,[105][106][107]: with, As Eq. 21 suggests, this new physics contribution to the T -parameter vanishes in the limit Since we are working in the exact alignment limit and a SFOPT demands m H ± = m A [64], hence T -parameter puts no bound on the scalar sector in our set-up.
The presence of the new VLFs shall also contribute to the T -parameter [108]: where The bound onŜ comes from a global fit: 10 3Ŝ = 0.0 ± 1.3 [109]. For S-parameter, we consider contribution only due to the VLFs as given by [37,42,108]: where g 2 is the SU (2) L gauge coupling. The expression for vacuum polarization for identical masses (at q 2 = 0) [37,42,108]: For two different masses (m i = m j ) the expression for vacuum polarization reads [37,42,108]: Note that all the divergences appearing in Eq. 25

Collider bounds
The LEP experiments have performed direct searches for charged Higgs. A combination of LEP data from searches in the τ ν and cs final states put a limit of m H ± > ∼ 80 Gev [110,111] [113][114][115][116]. Values of the charged Higgs masses for which such decay is kinematically allowed are excluded for tan β < ∼ 10 [40] in Type-I 2HDM. Here we would like to mention that in Type-I 2HDM tan β is unconstrained from Higgs signal strength in the strict alignment limit, that otherwise puts a strong limit [40,[117][118][119]. Flavour physics observable provide very strong constraints on charged Higgs mass. Inclusive b → sγ and more general b → s transitions lead to a robust exclusion of m H ± < 570 GeV [120] for Type-II 2HDM, while for Type-I it is excluded for tan β > ∼ 2 [111]. For Type-I 2HDM constraint from meson decay is rather weak, allowing Finally, We would like to highlight that LEP has set a lower limit on pair-produced charged heavy vector-like leptons: m ψ > ∼ 101.2 GeV at 95 % C.L. for ψ ± → νW ± final states [123]. For a SU (2) L singlet charged vector-like lepton the CMS search does not improve on the LEP limits. The limits for a heavy lepton doublet decaying to ∈ {e, µ} flavours are m L > ∼ 450 GeV [124]. In the case of decays to the τ flavour the limits are less stringent: m L > ∼ 270 GeV [124]. However, in our case the charged VLF ψ ± has dominant decay to the DM ψ 1 . As a result, the limits are less stringent and we only follow the LEP lmit in choosing our benchmark points for all the analyses.

Invisible decay constraints
When the DM mass m ψ 1 < m h 1 /2 or m ψ 1 < m Z /2, they can decay to a pair of the VLF DM (ψ 1 ). Higgs and Z invisible decays are precisely measured at the LHC [96]. Our model thus can be constrained from these measurements in the low mass range of the DM. Both the Higgs and Z invisible decays to DM are proportional to VLF mixing angle sin θ. In Appendix. A 1 we have computed the invisible decay width of Higgs and Z boson.

IV. DARK MATTER PHENOMENOLOGY
In this section we would like to elaborate on the DM phenomenology, where we show in detail the parameter space satisfying the PLANCK observed relic density by scanning over the free parameters of the model. Then we investigate how much of the relic density allowed parameter space also satisfies current direct detection bounds e.g, from XENON1T [11,12]. Finally, from the resulting parameter space satisfying relic abundance, direct search and existing theoretical and experimental bounds discussed earlier, we choose a few benchmark points for further analysis. All the relevant Feynman diagrams that contribute to the DM freeze-out are listed in Appendix. A 2.

A. Relic abundance of the dark matter
As we have already mentioned earlier, ψ 1 is the lightest VLF physical eigenstate which is odd under Z 2 , and hence a potential DM candidate in this model. The relic abundance of ψ 1 is mainly governed by the DM number changing annihilation and co-annihilation processes mediated by the where with n = n ψ 1 + n ψ 2 + n ψ ± and H is the Hubble parameter. In above equation, g ef f is defined as effective degrees of freedom, given by: where g 1 , g 2 and g 3 are the degrees of freedom of ψ 1 , ψ 2 and ψ ± respectively and We have implemented the model in LanHEP-3.3.2 [126] and the model files are then fed into micrOMEGAs-4.3.5 [127] for determining the relic abundance and direct detection cross-section for the DM. Before delving into the detailed parameter scan, we first start by looking into the variation of relic density with DM mass m ψ 1 for some fixed choices of two of the other free parameters: {∆m, sin θ}. Here we would like to clarify that for the entire analysis we have kept the masses of the new scalars fixed at: We perform a scan over a range of the parameter space: As discussed earlier in Sec. III, this choice of the scalar masses is motivated from the requirement of a SFOPT. Also, we would like to remind once more that we are strictly following the alignment limit, and hence the lightest CP-even scalar resembles the 125 GeV observed Higgs.
In is understandable, as larger sin θ gives rise to larger (co-)annihilation due to sin 2 θ dependence at the vertex for gauge-mediated processes and sin 2θ dependence for scalar mediated processes (due to proportionality to the Yukawa Y ). As a result, the relic abundance naturally decreases.
On the RHS in the bottom panel of Fig. 2 we show the same plot as that of the LHS but for larger where tan β = 1.3. Here we see, large ∆m is achieved for smaller sin θ. This can intuitively understood in the following way: the scalar mediated annihilation channels are essentially proportional to the Yukawa Y , which is proportional to both sin θ and ∆m. Hence a smaller sin θ requires a larger ∆m to produce the correct abundance. On the top right panel of Fig. 3 the same parameter space is shown for tan β = 5. We note the same pattern here except for the fact that 0.01 ≤ sin θ ≤ 0.1 region is more populated. This is a direct consequence of Eq. 17, which shows sin θ needs to be reduced as tan β increases (i.e., v 2 increases) to adjust the Yukawa Y such that the relic abundance is satisfied.
In order to understand the behaviour more intricately we have obtained a parameter space for a fixed sin θ in m ψ 1 -∆m plane as shown in the bottom right panel of Fig. 3. Here we have shown three different regions corresponding to under abundance (green), over abundance (red) and right relic (blue). As we move from left to right in this plot, we first encounter over abundant regions for small DM mass < ∼ 20 GeV. This is due to the lack of annihilation channels present for the DM to produce the right relic, as the only annihilation channels are to the light quarks. As we reach m ψ 1 ∼ 30 GeV, co-annihilation starts playing and as a consequence, the DM becomes under abundant. Still right relic abundance is not obtained as Y is small due to small ∆m and hence all the scalar mediated annihilations do not contribute significantly. Right relic is first obtained at 2 , due to the SM Higgs resonance. For DM mass ∼ 100 GeV as we move from lower ∆m to higher ∆m (from bottom to top), the DM is at first under abundant due to co-annihilation domination. Then right relic abundance is achieved as the co-annihilation is correctly tuned. Immediately after that there is an over abundant region for larger ∆m as co-annihilation loses its goodness and hence the effective annihilation cross-section (Eq. 28) becomes small. Note that, for a fixed DM mass > ∼ 100 GeV right relic abundance is reached twice: (a) Once for small ∆m, where right co-annihilation gives rise to observed relic and (b) for large ∆m, where Y is large enough to produce correct relic via scalar-mediated channels (to gauge-boson dominated final states) as shown in the bottom left panel of Fig. 3. Beyond m ψ 1 ∼ 500 GeV the parameter space is largely over abundant as the suppression due to 1/m 2 ψ 1 becomes significant, thus over producing the DM. Note that, right relic abundance is also obtained at the second Higgs resonance at m ψ 1 ∼ 150 GeV.

B. Direct detection of the dark matter
In the present framework, the DM exhibits spin-independent interactions with the nuclei induced at the tree-level by the mediation of CP-even states h 1 (SM Higgs-like) and h 2 (heavier Higgs), and also by the SM Z-boson exchange (as in Fig. 5). The relevant cross-section per nucleon reads [40]: where A is the mass number of the target nucleus, µ r = m ψ 1 m N m ψ 1 +m N is the DM-nucleus reduced mass and M is the spin-averaged DM-nucleus scattering amplitude given by: The effective couplings in Eq. 33 can be expressed as: with where f p,n T are nucleon form factors. For Z-mediated spin-independent direct detection, on the other hand, one can write the scattering cross-section per nucleon as [40]: with where f p,n are again suitable nucleon form factors 3 . For simplicity we can assume conservation of isospin i.e, f p /f n = 1. Here one should note that the Z-mediation also gives rise to spin-dependent direct search cross-section. But the order of magnitude of such spin-dependent cross-section being extremely small compared to that of spin-independent ones we choose to ignore that. 4 .
In the top left panel of Fig. 6 we have illustrated the relic density allowed parameter space that survives present spin-independent direct detection bound from XENON1T for tan β = 1.3. We see a moderate range of sin θ's are allowed by direct search: 0.01 < ∼ sin θ < ∼ 0.3. This is expected as the direct detection cross-section is proportional to sin 2 θ for scalar-mediation and sin 4 θ for gaugemediation. As a consequence, smaller mixing should give rise to smaller σ SI , making the DM parameter space more viable from direct search bound. The presence of the second Higgs helps in keeping the VLF mixing within a moderate limit, unlike the case in [32] where the bound on the mixing is even more stringent due to the presence of only one (SM) Higgs. This is a consequence of sin 2 α m 2 H suppression due to the heavier Higgs and small scalar mixing in case of scalar-mediated elastic scattering. This is also possible due to some cancellation between the Higgs-mediated diagrams leading to a destructive interference that allows one to choose sin θ as large as ∼ 0.3 without getting disallowed by the direct search exclusion. Also note here, most of the allowed parameter space lies just above the neutrino floor [128] and hence can still be probed by the future direct search experiments with improved sensitivity. In summary, constraints from the requirement of right relic abundance, together with the direct search exclusion limit allows the VLF mixing to vary within a range of 0.01 < ∼ sin θ < ∼ 0.3 for a DM mass starting from around 100 GeV upto 3 TeV.
The presence of the second Higgs helps the model to evade present direct search bound and allows the parameter space to fit just above the neutrino floor, leaving the window open to either get discovered or get discarded from the very next limit on spin-independent direct search. Top right panel of Fig. 6 shows the same with tan β = 5.
In the bottom left panel of Fig. 6 we show the residual parameter space satisfying both relic abundance and direct detection bounds in m ψ 1 -∆m plane with respect to the variation of the VLF mixing sin θ for tan β = 1.3. This clearly shows that ∆m < ∼ 1.5 TeV (for DM mass ∼ 500 GeV) in order to abide by both relic abundance and spin-independent direct detection bounds for sin θ < ∼ 0.3.
For larger tan β, shown in the right panel of Fig. 6, the bound on ∆m is bit more relaxed, which allows it to ∼ 2.5 TeV but for a larger DM mass. However, the VLF mixing is rather restricted and can be as large as sin θ ∼ 0.1, which helps in suitably choosing the Yukawa Y via Eq. 17. The bound on ∆m is crucial as larger ∆m results in larger missing energy, which, in turn helps the model to be separated from the SM background at the colliders as we shall explain in Sec. V. For smaller ∆m the model can still be found at the colliders via stable charged track signature. Before moving on to the collider analysis we have tabulated some of the benchmark points (BP) in Tab. II. These BPs satisfy bounds from relic abundance, direct detection and also those arising from EWPO. The BPs are listed in the decreasing order of ∆m, where BP6 has the minimum ∆m.
We have also kept the DM mass < 2 TeV to have a sizeable ψ ± production cross-section. As we shall show in Sec. V, BP(1-5) can be distinguished at the LHC as they produce huge missing energy due to large ∆m that helps to separate them from the SM background. On the other hand, due to very small ∆m (specifically as ∆m < m W ), BP6 can only produce a displaced vertex at the collider or can be searched at the ILC via missing energy excess as shown in [36,37,41].

V. COLLIDER PHENOMENOLOGY
The detailed study of collider signature for vector-like fermions can be found in [32,36,37,41]. As we have already seen, due to the presence of the second Higgs doublet large sin θ can be achieved satisfying both relic density and direct detection. As a result one needs not be confined in small ∆m and a large ∆m is also viable. Such large ∆m's are favorable in order to distinguish this model at the collider from the SM background [41,42]. It is to be noted that the charged component of SU (2) L doublet VLF can be produced at the LHC via SM Z and photon mediation. The charged VLF can further decay via on-shell and/or off-shell W (depending on whether ∆M > ∼ 80 GeV or ∆m < ∼ 80 GeV) to the following final states: • Hadronically quiet opposite sign dilepton (OSD) with missing energy ( + − + / E T ).
• Four jets plus missing energy (jjjj + / E T ). We shall focus only on the leptonic final states (hadronically quiet dilepton) as they are much cleaner compared to others. The Feynman graph for such a process is depicted in the LHS of Fig. 7.
In the RHS of Fig. 7 we have shown the variation of pair production cross-section of the charged component of the VLF with VLF mass at √ s = 14 TeV. As one can see, the production cross-section decreases with increase in the charged VLF mass showing the usual nature. We have also shown the position of different BPs on the same plot. As one can notice, BP6 has the highest production cross-section while BP2 has the least. This is evident from the fact that for BP6 ∆m = 10 GeV, giving rise to charged VLF mass of m ψ ± = 276 GeV. On the other hand BP2 has a larger ∆m but highest DM mass, which gives rise to m ψ ± = 1803 GeV much larger than that for BP4.

A. Object reconstruction and simulation details
As already mentioned, we implemented this model in LanHEP-3.3.2 and the parton level events are generated in CalcHEP-3.7.3 [129]. Those events are then fed to PYTHIA-6.4 [130] for showering and hadronization. The dominant SM backgrounds that can imitate our final state are generated in MADGRAPH-2.6.6 [131], and the corresponding production cross-sections are multiplied with appropriate K-factor [132] in order to match with the next to leading order (NLO) cross-sections.
For all cases we have used CTEQ6l as the parton distribution function (PDF) [133].
where the sum runs over all visible objects that include the leptons, jets and the unclustered components.

v) Invariant dilepton mass (m ):
We can construct the invariant dilepton mass variable for two opposite sign leptons by defining: Invariant mass of OSD events, if created from a single parent, peak at the parent mass, for example, Z boson. As the signal events (Fig. 7) do not arise from a single parent particle, invariant mass cut plays key role in eliminating the Z mediated SM background.
vi) H T : H T is defined as the scalar sum of all isolated jets and lepton p T 's: For our signal the sum only includes the two leptons that are present in the final state.
We shall use different cuts on these observables depending on their distribution patterns to separate the signal from the SM backgrounds. Thus we can predict the significance as a function of the integrated luminosity. These are discussed in the following sections.

B. Event rates and signal significance
In Fig. 8   for larger ∆m the signal distributions are well separated from that of the background. This is due to the fact that the peak of the MET distribution is determined by how much of p T is being carried away by the missing particle (i.e, the DM), which in turn depends on the mass difference of charged and neutral component of the VLF i.e, ∆m. Hence for larger ∆m the DM carries away most of the p T making the distribution much flatter, while for smaller ∆m the distribution peaks up at lower value as the produced DM particles are not boosted enough. As a consequence, BP1 and BP2 have the most flattened distribution, while BP (3,4,5) are increasingly overwhelmed by the SM background.
Since BP6 has the smallest ∆m, we refrain from showing this in the plots as BP6 will be inseparable from the SM backgrounds. From these distributions it is quite evident that with a cut on MET This can be translated into the signal significance for this model which is shown in Fig. 9  even with a luminosity as high as 3000 fb −1 . A large ∆m in one way helps to distinguish BP (2,4) from the SM background, but due to very small cross-section at the production level the final state cross-section becomes even smaller. This makes BP(2,4) least significant. For BP(1,3,5), on the other hand, ∆m is optimum such that they can be separated from the background because of their larger missing energy, and at the same time the cross-section is large enough so that a 5σ significance can be achieved. Thus, we see, BP3 and BP5 reach a discovery limit at a lower luminosity ∼ 500 fb −1 , while BP1 needs ∼ 1000 fb −1 for a 5σ reach. Finally, we would like to mention that for all benchmarks with ∆m < ∼ m W (e.g, BP6) one can find a stable charged track due to the off-shell decay of the heavy charged VLF ψ ± via W -boson or they can also be probed at the ILC with a lower cut on MET. These scenarios have already been thoroughly investigated in [36,41], hence we do not discuss here further.

WAVE SIGNALS
In this section we would like to show the possibility of generation of stochastic gravitational wave (GW) from a strong first-order phase transition (SFOPT). The frequency of such GWs are well within the reach of the proposed GW detectors. The occurrence of a SFOPT and subsequent GW generation in context with 2HDM has already been thoroughly studied [64,65]. In the context of our present model we explore how such a detectable GW signal can be an alternate search strategy for singlet-doublet DM. The dynamics of the SFOPT is completely determined by the parameters of the scalar potential as we shall see in the following sections. It is interesting to note that the choice of the scalar potential parameters agrees well with both the DM phenomenology and the GW generation, thus giving us a handle to probe the dark sector beyond DM and collider search experiments.

A. Finite Temperature Effective Potential
In order to explore the electroweak phase transition (EWPT) in 2HDM we need to include temperature corrections with the tree-level potential. In general, the finite temperature effective potential at a temperature T can be expressed as [134] where V tree , V T =0 1−loop and V T =0 1−loop are the tree-level potential at zero temperature, the Coleman-Weinberg one-loop effective potential at zero temperature and one-loop effective potential at finite temperature respectively. The tree-level potential V tree can be obtained from Eq. 5 by replacing the fields Φ 1 and Φ 2 with their classical fields v 1 and v 2 which is given by The Coleman-Weinberg one-loop effective potential at zero temperature V T =0 1−loop can be written as [134,135] where y t , g and g are the top Yukawa coupling, SU (2) L and U (1) Y gauge couplings of the SM respectively.
The field-dependent squared masses m 2 i at T = 0 for the scalar bosons can be obtained by diagonalizing the following matrices Here we have applied Landau gauge, where the Goldstones are massless at zero temperature (T = 0) but at finite temperature (T = 0) they acquire a mass [59].
where the function J ± are In the finite temperature effective potential V T =0 1−loop we include the temperature corrected terms to the boson masses by following Daisy resummation method [136]. In the Daisy resummation method the thermal masses are [137], [138][139][140]

B. Gravitational Wave from SFOPT
The central idea of first order phase transition (FOPT) is the bubble nucleation of a true vacuum state (from several metastable states) at a temperature commonly known as the nucleation temperature. The bubbles produced in this process can be of different sizes: small and large. The smaller bubbles tend to collapse, whereas the larger bubbles tend to expand after attaining the criticality. These bubbles of critical size then collide with each other and their spherical symmetry is thus broken. This initiates the phase transition and subsequent production of the GW. The bubble nucleation rate per unit volume at a temperature T can be expressed as [141] where Γ 0 (T ) ∝ T 4 and S 3 (T ) denotes the Euclidean action of the critical bubble [141]: where V ef f is the effective finite temperature potential (Eq. (42)). Bubble nucleation occurs at the nucleation temperature T n which satisfies the condition: S 3 (T n ) /T n ≈ 140 [134].
As mentioned in the Sec. I, GWs are produced from the FOPT majorly via three mechanisms namely, bubble collisions [66][67][68][69][70][71][72], sound wave [73][74][75][76] and turbulence in the plasma [77][78][79][80][81]. The total GW intensity Ω GW h 2 as a function of frequency can be expressed as the sum of the contributions from the individual components [66][67][68][69][70][71][72][73][74][75][76][77][78][79][80][81]: The component from the bubbles collision Ω col h 2 is given by (for an analytic and more accurate derivation see [142,143]): where the parameter β β = HT d dT where T n is the nucleation temperature and H n is the Hubble parameter at T n . The most general expression of the bubble wall velocity v w can be written as 5 [71,145] v w = 1/ √ 3 + α 2 + 2α /3 The parameter κ in Eq. (58) is the fraction of latent heat deposited in a thin shell which can be expressed as: with [85,92] α ∞ = 30 where v n represents the vacuum expectation value of Higgs at T n and m W , m Z and m t are the masses of W, Z and top quarks respectively. The parameter α , which is defined as the ratio of vacuum energy density ρ vac released by the electroweak phase transition to the background energy density of the plasma ρ rad * at T n , has the form: with and The quantity f col in Eq. (58) is the peak frequency produced by the bubble collisions and reads: The sound wave (SW) component of the gravitational wave (Eq. (57)) is given by where κ v is the faction of latent heat transformed into the bulk motion of the fluid which can be expressed as In Eq. (67) f SW denotes the peak frequency produced by the sound wave mechanisms which has the following form To check the contribution of sound wave component to the total GW intensity we need to estimate the suppression factor HR * /Ū f , whereŪ f denotes the root-mean-square (RMS) fluid velocity and R * denotes the mean bubble separation [85,146,147]. If the calculated suppression factor HR * /Ū f of a given model comes out to be > 1 , then the sound wave lasts more than a Hubble time, otherwise it is an overestimate to the GW signal.  Finally, the component from the turbulence in the plasma Ω turb h 2 is given by: where = 0.1 and f turb denotes the peak frequency produced by the turbulence mechanism which takes the form In Eq. (70) the parameter h * can be written as h * = 16.5 × 10 −6 Hz T n 100 GeV g * 100 Eq. (57)-(72) are used for calculating the gravitational wave intensity. In order to study the phase transition properties and production of GW in the present DM model, we choose four BPs (Table   V) from the viable model parameter space. An EWPT takes place at the nucleation temperature where a high phase and a low phase is separated by a potential barrier. The strength of the phase transition depends on the ratio of the VEVs v c = √ < v 1 > 2 + < v 2 > 2 , measured at the critical temperature, to the critical temperature T c at which two degenerate minima exist. If the condition ξ = v c /T c > 1 is satisfied then the strong first-order phase transition (SFOPT) is said to be occured [64]. We calculate the GW intensity using Eq. (57)- (72). For computing the intensity we need to first estimate the thermal parameters related to FOPT. In order to find them we have used Cosmotransition package [134]. The tree level potential (Eq. 43) is given as an input to this package, and resulting thermal parameters are tabulated in Table VI. The GW intensity mainly depends on several factors e.g, nucleation temperature T n , bubble wall velocity v w , strength of the FOPT α and the parameter β . As mentioned in Sec. VI B, the sound wave contribution to the total GW intensity depends on the suppression factor HR * /Ū f depending on whether it lasts more than a Hubble time or not. We estimate the suppression factor HR * /Ū f following [85,146,147] and found it to be < 1 for all the BPs. Due to this fact, following [85,146,147], we include the suppression factor HR * /Ū f to the sound wave component of the GW intensity. In Figure 10 we have plotted and compared the GW intensities for the chosen benchmarks (Table V)  DECIGO, aLIGO, aLIGO+ and LISA following Ref. [149,150]. The frequencies at which the GW intensities acquire maximum value are 10 −2 Hz, 7 × 10 −3 Hz, 3.80 × 10 −2 Hz and 3.40 × 10 −2 Hz for BPI, BPII, BPIII, and BPIV respectively. As it is evident from Figure 10, the GW intensities for all the BPs (BPI-BPIV) lie within the sensitivity curves of ALIA, BBO and DECIGO. The upshot is thus to note that fact that the benchmark values of m 12 and tan β for DM phenomenology agrees well with that of a successful SFPOT leading to production of detectable GW signal.

VII. CONCLUSIONS
In this paper we have proposed a singlet-doublet fermionic dark matter (DM) model, where the lightest fermion (odd under an imposed discrete symmetry Z 2 ), emerging as a singlet-doublet admixture of vectorlike fermions (VLF), can be a viable DM candidate (ψ 1 ). We extend the model with a second Higgs doublet where the second Higgs is odd under another discrete symmetry Z 2 . The imposition of two different discrete symmetries is necessary to prevent the occurrence of the flavour changing neutral current (FCNC) at tree-level, while allowing a strong first order phase transition within a consistent DM framework. The DM in this case is a weakly interacting massive particle that undergoes freeze-out to yield the PLANCK observed relic abundance. This is achieved via annihilation of the DM with itself, also via its co-annihilation with its massive component (ψ 2 ) and with the charged component (ψ ± ). As the doublet carries a SU (2) L charge, hence on top of scalar mediation, the annihilation channels are also SM gauge mediated. The presence of the Z-mediated direct search puts a strong bound on the model parameter space allowing the VLF mixing sin θ < ∼ 0.3 both for tan β = 1. 3 and tan β = 5 for DM mass up to > ∼ 1 TeV. The presence the second Higgs makes the direct detection bounds less stringent. This is typically due to (a) sin 2 α/m 2 H suppression from the heavier Higgs with small scalar mixing and (b) some destructive interference between the two scalar mediated diagrams (the so called "blind spot") that offers some breathing space in the direct detection parameter space. Thus one can still achieve a moderate sin θ in contrast with singlet-doublet models with only SM Higgs. The model thus lives over a large parameter space satisfying both relic abundance and spin-independent direct search bounds.
We have then explored possible signatures that this model can give rise to at the LHC. As the leptonic channels are cleaner, we have studied the hadronically quiet dilepton final states (HQ2L) where we see, a substantial signal significance is achievable (for an integrated luminosity L ∼ 300 fb −1 ) by a judicious choice of cuts on the missing energy (MET) and H T . This is again possible because of the presence of the second Higgs doublet that allows large sin θ satisfying both relic abundance and direct search, which, in turn, allows a large ∆m ≡ m ψ ± − m ψ 1 ∼ 1 TeV. Larger ∆m is the key to distinguish the signal from the dominant SM backgrounds exploiting hard cuts on MET and H T . For ∆m < ∼ m W the model may be probed via displaced vertex signature due to the off-shell decay of the charged VLF to SM leptons and neutrino.
We finally have looked into the possibility of getting gravitational wave (GW) via a strong firstorder phase transition (SFOPT) due to the extended scalar sector. We have found for both tan β = {1.3, 5} the model is capable of producing detectable GW signal modulo we tune the parameter m 12 accordingly. The probability of getting a SFOPT is thus very sensitive to the choice of m 12 . We see, for some benchmark values of m 12  Note that for none of the benchmark points in Tab. II, except BP3, either of the constraints from Higgs invisible decay branching or Z-boson invisible decay branching is applicable. Since for BP3 the DM mass is 60 GeV, it is possible for the SM Higgs to decay to a pair of ψ 1 . However, due to small VLF mixing such a decay is well within the measured invisible decay rate of SM Higgs. 11. Annihilation (i = j) and Co-annihilation (i = j) type number changing processes for Vector like fermionic DM in the model. Here i, j, k = 1, 2; a, b, c = 1, 2 and f stands for SM fermions.