Dark matter, electroweak phase transition and gravitational wave in the type-II two-Higgs-doublet model with a singlet scalar field

In the framework of type-II two-Higgs-doublet model with a singlet scalar dark matter $S$, we study the dark matter observables, the electroweak phase transition, and the gravitational wave signals by such strongly first order phase transition after imposing the constraints of the LHC Higgs data. We take the heavy CP-even Higgs $H$ as the only portal between the dark matter and SM sectors, and find the LHC Higgs data and dark matter observables require $m_S$ and $m_H$ to be larger than 130 GeV and 360 GeV for $m_A=600$ GeV in the case of the 125 GeV Higgs with the SM-like coupling. Next, we carve out some parameter space where a strongly first order electroweak phase transition can be achieved, and find benchmark points for which the amplitudes of gravitational wave spectra reach the sensitivities of the future gravitational wave detectors.


I. INTRODUCTION
The weakly interacting massive particle is a primary candidate for dark matter (DM) in the present Universe. Many extensions of SM have been proposed to provide a candidate of DM, and one simple extension is to add a singlet scalar DM to the type-II two-Higgs-doublet model (2HDM) [1][2][3]. The type-II 2HDM (2HDMIID) contains two CP-even states, h and H, one neutral pseudoscalar A, two charged scalars H ± , and one CP-even singlet scalar S as the candidate of DM [4][5][6][7][8][9][10][11][12][13][14].
In the type-II 2HDM model, the Yukawa couplings of the down-type quark and lepton can be both enhanced by a factor of tan β. Therefore, the flavor observables and the LHC searches for Higgs can impose strong constraints on type-II 2HDM model. In the 2HDMIID, the two CP-even states h and H may be portals between the DM and SM sectors, and there are plentiful parameter space satisfying the direct and indirect experimental constraints of DM. The scalar potential of 2HDMIID contains the original one of type-II 2HDM and one including DM. For appropriate Higgs mass spectrum and coupling constants, the type-II 2HDM can trigger a strong first-order electroweak phase transition (SFOEWPT) in the early universe [15][16][17][18], which is required by a successful explanation of the observed baryon asymmetry of the universe (BAU) [19] and can produce primordial gravitational wave (GW) signals [20].
In this paper, we first examine the parameter space of the 2HDMIID using the recent LHC Higgs data and DM observables. After imposing various theroretial and experimental constraints, we analyze whether a SFOEWPT is achievable in the 2HDMIID, and discuss the resultant GW signals and its detectability at the future GW detectors, such as LISA [21], Taiji [22], TianQin [23], Big Bang Observer (BBO) [24], DECi-hertz Interferometer GW Observatory (DECIGO) [24] and Ultimate-DECIGO [25].
Our work is organized as follows. In Sec. II we will give a brief introduction on the 2HDMIID. In Sec. III and Sec. IV, we show the allowed parameter space after imposing the limits of the LHC Higgs data and DM observables. In Sec. V, we examine the parameter space leading to a SFOEWPT and the corresponding GW signal. Finally, we give our conclusion in Sec. VI.

II. TYPE-II TWO-HIGGS-DOUBLET MODEL WITH A SCALAR DARK MAT-TER
The scalar potential of 2HDMIID is given as [26] V tree = m 2 Here we discuss the CP-conserving model in which all λ i , κ i and m 2 12 are real. The S is a real singlet scalar field, and Φ 1 and Φ 2 are complex Higgs doublets with hypercharge Y = 1: Where v 1 and v 2 are the electroweak vacuum expectation values (VEVs) with v 2 = v 2 1 + v 2 2 = (246 GeV) 2 , and the ratio of the two VEVs is defined as tan β = v 2 /v 1 . The linear and cubic terms of the S field are forbidden by a Z 2 symmetry, under which S → −S. The S is a possible DM candidate since it does not acquire a VEV. After spontaneous electroweak symmetry breaking, the remaining physical states are three neutral CP-even states h, H, and S, one neutral pseudoscalar A, and two charged scalars H ± .
The Yukawa couplings of the neutral Higgs bosons normalized to the SM are given by The charged Higgs has the following Yukawa interactions, where i, j = 1, 2, 3.
The neutral Higgs boson couplings with the gauge bosons normalized to the SM are given by where V denotes Z or W . In the type-II 2HDM, the 125 GeV Higgs is allowed to have the SM-like coupling and wrong sign Yukawa coupling, Therefore, in this paper we take the heavy CP-even Higgs H as the only portal between the DM and SM sectors, and focus on the case of the 125 GeV with the SM-like couping.
The S, T , and U oblique parameters give the stringent constraints on the mass spectrum of Higgses of type-II 2HDM [28][29][30]. One of m A and m H is around 600 GeV, and another is allowed to have a wide mass range including low mass [29,30]. Therefore, we fix m A = 600 GeV to make the portal H to have a wide mass range.
In our calculation, we consider the following observables and constraints: (1) Theoretical constraints. The scalar potential of the model contains one of the type-II 2HDM and one of the DM sector. The vacuum stability, perturbativity, and tree-level unitarity impose constraints on the relevant parameters, which are discussed in detail in Refs. [8,9]. Here we employ the formulas in [8,9] to implement the theoretical constraints. Compared to Refs. [8,9], there are additional factors of 1 2 in the κ 1 term and the κ 2 term of this paper. In addition, we require that the potential has a global minimum at the point of (< h 1 >= v 1 , < h 2 >= v 2 , < S 1 >= 0).
(2) The oblique parameters. The S, T , U parameters can impose stringent constraints on the mass spectrum of Higgses of 2HDM. We use 2HDMC [31] to calculate the S, T , U parameters. Taking the recent fit results of Ref. [32], we use the following values of S, T , U , S = 0.02 ± 0.10, T = 0.07 ± 0.12, U = 0.00 ± 0.09.
The correlation coefficients are (3) The flavor observables and R b . We employ SuperIso-3.4 [33] to calculate Br(B → X s γ), and ∆m Bs is calculated following the formulas in [34]. Besides, we include the constraints of bottom quarks produced in Z decays, R b , which is calculated following the formulas in [35,36].
These samples correspond to be within the 2σ range in any two-dimension plane of the model parameters when explaining the Higgs data.
(5) The exclusion limits of searches for additional Higgs bosons. We use HiggsBounds-4.3.1 [38,39] to implement the exclusion constraints from the neutral and charged Higgs searches at LEP at 95% confidence level.
Because the b-quark loop and top quark loop have destructive interference contributions to gg → A production in the type-II 2HDM, the cross section decreases with an increase of tan β, reaches the minimum value for the moderate tan β, and is dominated by the b-quark loop for enough large tan β. In addition to tan β and m H , the cross section of gg → H depends on sin(β − α). We employ SusHi to compute the cross sections for H and A in the gluon fusion and bb-associated production at NNLO in QCD [40]. The cross sections of H via vector boson fusion process are deduced from results of the LHC Higgs Cross Section Working Group [41]. We employ 2HDMC to calculate the branching ratios of the various decay modes of H and A. The searches for the additional Higgs considered by us are listed in Tables I and II. The LHC searches for H ± can not impose any constraints on the model for m H ± > 500 GeV and 1 ≤ tan β ≤ 25 [42]. Therefore, we do not consider the constraints from the searches for the heavy charged Higgs.

B. Results and discussions
In Fig as only portal between DM and SM sectors, and the decay H → SS opens for 2m S < m H .
The decay mode possibly affects the allowed parameter space, but the constraints of the DM observables have to be simultaneously considered. Here we temporarily assume 2m S > m H , and close the H → SS decay mode. In the next section, the effects of H → SS will be considered by combining the DM observables.
In increases with m H .

IV. THE DARK MATTER OBSERVABLES
We use FeynRules [87] to generate the model file, which is called by micrOMEGAs [88] to calculate the relic density. In our scenario, the elastic scattering of S on a nucleon receives the contributions of the process with t-channel exchange of H, and the spin-independent cross section between DM and nucleons is given by [89] σ p(n) = µ 2 where with

V. ELECTROWEAK PHASE TRANSITION AND GRAVITATIONAL WAVE
The phase transition can proceed in basically two different ways. In a first-order phase transition, at the critical temperature T C , the two degenerate minima will be at different points in field space, typically with a potential barrier in between. For a second order (crossover) transition, the broken and symmetric minimum are not degenerate until they are at the same point in field space. In this paper we focus on the SFOEWPT, which is required by a successful explanation of the observed BAU and can produce primordial GW signals.
The masses of gauge boson are given We neglect the contributions of light fermions, and only consider the masses of top quark and bottom quark, where v . Now we study the effective potential with thermal correction. The thermal effective potential V ef f in terms of the classical fields (h 1 , h 2 , S 1 ) is composed of four parts: Where V 0 is the tree-level potential, V CW is the Coleman-Weinberg potential, V CT is the counter term, V T is the thermal correction, and V ring is the resummed daisy corrections. In this paper, we calculate V ef f in the Landau gauge.
We obtain the tree-level potential V 0 in terms of their classical fields (h 1 , h 2 , S 1 ) The Coleman-Weinberg potential in the MS scheme at 1-loop level has the form [93]: where i = h, H, A, H ± , S, G, G ± , W ± , Z, t, b, and s i is the spin of particle i. Q is a renormalization scale, and we take Q 2 = v 2 . The constants C i = 3 2 for scalars or fermions and C i = 6 5 for gauge bosons. n i is the number of degree of freedom, n h = n H = n G = n A = 1, n H ± = n G ± = 2, n W ± = 6, n Z = 3, With V CW being included in the potential, the minimization conditions of scalar potential in Eq. (19) and the CP-even mass matrix will be shifted slightly. To maintain the minimization conditions at T=0, we add the so-called "counter-terms" where the relevant coefficients are determined by which are evaluated at the EW minimum of {h 1 = vc β , h 2 = vs β , S 1 = 0} on both sides. As a result, the VEVs of h 1 , h 2 , S 1 and the CP-even mass matrix will not be shifted.
It is a well-known problem that the second derivative of the Coleman-Weinberg potential at T = 0 suffers from logarithmic divergences originating from the vanishing Goldstone masses. To solve the divergence problem, we take a straightforward approach of imposing an IR cut-off at m 2 IR = m 2 h for the masses of Goldstone boson of the divergent terms, which gives a good approximation to the exact procedure of on-shell renormalization, as argued in [16].
The thermal contributions V T to the potential can be written as [94] where i = h, H, A, H ± , S, G, G ± , W ± , Z, t, b, and the functions J B,F are Finally, the thermal corrections with resumed ring diagrams are given [95,96] where i = h, H, A, H ± , S, G, G ± , W ± L , Z L , γ L . The W ± L , Z L , and γ L are the longitudinal gauge bosons with n W ± L = 2, n Z L = n γ L = 1. The thermal Debye massesM 2 i (h 1 , h 2 , S 1 , T ) are the eigenvalues of the full mass matrix, where X = P, A, C. Π X are given by The physical mass of the longitudinally polarized W boson is The physical mass of the longitudinally polarized Z and γ boson with

B. Calculation of electroweak phase transition and gravitational wave
In a first-order cosmological phase transition, bubbles nucleate and expand, converting the high-temperature phase into the low-temperature one. The bubble nucleation rate per unit volume at finite temperature is given by [97][98][99] where A(T ) ∼ T 4 is a prefactor and S E is the Euclidean action At the nucleation temperature T n , the thermal tunneling probability for bubble nucleation per horizon volume and per horizon time is of order one, and the conventional condition is There are two key parameters characterizing the dynamics of the EWPT, β and α. β describes roughly the inverse time duration of the strong first order phase transition, where H n is the Hubble parameter at the bubble nucleation temperature T n . α is defined as the vacuum energy released from the phase transition normalized by the total radiation energy density ρ R at T n , where g * is the effective number of relativistic degrees of freedom. We use the numerical package CosmoTransitions [100] and PhaseTracer [101] to analyze the phase transition and computes quantities related to cosmological phase transition.
In a radiation-dominated Universe, there are three sources of GW production at a EWPT: bubble collisions, in which the localized energy density generates a quadrupole contribution to the stress-energy tensor, which in turn gives rise to GW, plus sound waves in the plasma and magnetohydrodynamic (MHD) turbulence. The total resultant energy density spectrum can be approximately given as, Recent studies show that the energy deposited in the bubble walls is negligible, despite the possibility that the bubble walls can run away in some circumstances [102]. Therefore, although a bubble wall can reach relativistic speed, its contribution to GW can generally be neglected [103,104]. Therefore, in the following discussions we do not include the contribution from bubble collision Ω col .
The GW spectrum from the the sound waves can be obtained by fitting to the result of numerical simulations [105], where f sw is the present peak frequency of the spectrum, Hz . (40) v w is the wall velocity, and the factor κ v is the fraction of latent heat transformed into the kinetic energy of the fluid. κ v and v w are difficult to compute, and involves certain assumptions regarding the dynamics of the bubble walls. On the other hand, successful electroweak baryogenesis scenarios favor lower wall velocity v w ≤ 0.15 − 0.3 [106], which allows the effective diffuse of particle asymmetries near the bubble wall front. In Ref. [107], however, it is pointed out that the relevant velocity for electroweak baryogenesis is not really v w , but the relative velocity between the bubble wall and the plasma in the deflagration front.
As a result, the electroweak baryogenesis is not necessarily impossible even in the case with large v w . Therefore, in this paper we take two different cases of v w and κ v [108]: • For small wall velocity: v w = 0.3 and • For very large wall velocity: v w = 0.9 and Considering Kolmogorov-type turbulence as proposed in Ref. [109], the GW spectrum from the MHD turbulence has the form [110,111], with the red-shifted Hubble rate at GW generation h n = 1.65 × 10 −5 T n 100GeV g * 100 1 6 Hz.
The peak frequency f turb is given by Hz .
The energy fraction transferred to the MHD turbulence κ turb can vary between 5% to 10% of κ v [105]. Here we take κ turb = 0.1κ v . For both sound wave and turbulence contribution as shown in Eq. (40) and Eq. (43), the amplitudes of the GW spectra are proportional to v w and the peak frequencies shift as 1/v w . Therefore, one changes in the wall velocity approximately have an order one effect on the spectrum and peak frequencies.

C. Results and discussions
The strength of the electroweak phase transition is quantified as with v c = √ < h 1 > 2 + < h 2 > 2 at critical temperature T c . The global minimum of potential has < A >= 0 because of the CP-conserving case. In order to avoid washing out the baryon number generated during the phase transition, a SFOEWPT is required and the conventional condition is ξ c ≥ 1.
After imposing the constraints of "pre-LHC", the LHC Higgs data, the relic density, XENON1T, and Fermi-LAT, we scan over the parameter space in the previous selected scenario. We find some surviving samples which can achieve a SFOEWPT, and these samples are projected in Fig. 4 and Fig. 5. For all the surviving samples, at T c the two degenerate minima of potential are respectively at (< h 1 >, < h 2 >, 0) and (0, 0, 0). In the process of EWPT, < S 1 > always has no VEV. an increase of tan β. There is a relative strong correlation between m 2 12 and tan β, and m 2 12 is imposed upper and lower bounds for a given tan β. With an increase of tan β, m 2 12 is stringently restricted by the theoretical constraints and the LHC Higgs data, leading that it is difficult to achieve a SFOEWPT. Thus, most of samples lie in the region of small tan β.
The requirement of a SFOEWPT is not sensitive to m S , and disfavors | λ H |> 0.3. Now we examine two key parameters α and β/H n which characterize the dynamics of the SFOEWPT, and govern the strength of GW spectra. A larger α and smaller β/H n can lead to stronger GW signals. In addition to the conditions of the successful bubble nucleations, we require with v n = √ < h 1 > 2 + < h 2 > 2 at the nucleation temperature T n . In fact, this is a more precise condition of SFOEWPT than ξ c ≥ 1. Also note that there generically exists a difficulty for solving bounce solution in a very thin-walled bubble, including the package CosmoTransitions [112]. Therefore, we will neglect the samples with very thin-walled bubble. Consider the constraints discussed above, we find some surviving samples, and the corresponding α and β/H n are shown in Fig. 6.
The β/H n may characterize the inverse time duration of the EWPT. A small β/H n means a long EWPT, and gives strong GW signals. For the GW coming from the sound waves in the plasma, the GW signal will continue being generated and the energy density of the GW  is thus proportional to the duration of the EWPT if the mean square fluid velocity of the plasma is non-negligible [105]. In addition, a large β/H n can enhance the peak frequency of the GW spectra. The parameter α describes the amount of energy released during the EWPT, and therefore a large α leads to strong GW signals.
We pick out two benchmark points (BPs), and examine the corresponding GW spectra.

VI. CONCLUSION
We examine the status of the 2HDMIID confronted with the recent LHC Higgs data, the DM observables and SFOEWPT, and discuss the detectability of GW at the future GW A SFOEWPT can be achieved in the many regions of m H < 500 GeV and m A = 600 GeV, favors a small tan β, and is not sensitive to the mass of DM. We find the benchmark points for which the predicted GW spectra can reach the sensitivities of LISA, TianQin, BBO, DECIGO, and UDECIGO. [54] ATLAS collaboration, "Search for a high-mass Higgs boson decaying to a pair of W bosons in pp collisions at √ s = 13 TeV with the ATLAS detector," ATLAS-CONF-2016-074.
[55] ATLAS Collaboration, "Search for diboson resonance production in the νqq final state using p p collisions at √ s = 13 TeV with the ATLAS detector at the LHC," ATLAS-CONF-2016-062.
[56] ATLAS Collaboration, "Search for WW/WZ resonance production in νqq final states in pp collisions at √ s = 13 TeV with the ATLAS detector," arXiv:1710.07235.
[58] CMS Collaboration, "Search for a heavy Higgs boson decaying to a pair of W bosons in proton-proton collisions at √ s = 13 TeV," arXiv:1912.01594. [62] ATLAS Collaboration, "Study of the Higgs boson properties and search for high-mass scalar resonances in the H → ZZ * → 4 decay channel at √ s = 13 TeV with the ATLAS detector," ATLAS-CONF-2016-079.
[64] ATLAS Collaboration, "Searches for heavy ZZ and ZW resonances in the qq and ννqq final states in pp collisions at √ s = 13 TeV with the ATLAS detector," arXiv:1708.09638. [68] ATLAS Collaboration, "Search for pair production of Higgs bosons in the bbbb final state using proton−proton collisions at √ s = 13 TeV with the ATLAS detector," ATLAS-CONF-2016-049.
[69] CMS Collaboration, "Search for a massive resonance decaying to a pair of Higgs bosons in the four b quark final state in proton-proton collisions at √ s = 13 TeV," arXiv:1710.04960.
[70] CMS Collaboration, "Search for Higgs boson pair production in events with two bottom quarks and two tau leptons in proton-proton collisions at √ s = 13 TeV," arXiv:1707.02909.
[73] ATLAS Collaboration, "Reconstruction and identification of boosted di-τ systems in a search for Higgs boson pairs using 13 TeV proton-proton collision data in ATLAS," arXiv:2007.14811.