Phase transition gravitational waves from pseudo-Nambu-Goldstone dark matter and two Higgs doublets

We investigate the potential stochastic gravitational waves from first-order electroweak phase transitions in a model with pseudo-Nambu-Goldstone dark matter and two Higgs doublets. The dark matter candidate can naturally evade direct detection bounds, and can achieve the observed relic abundance via the thermal mechanism. Three scalar fields in the model obtain vacuum expectation values, related to phase transitions at the early Universe. We search for the parameter points that can cause first-order phase transitions, taking into account the existed experimental constraints. The resulting gravitational wave spectra are further evaluated. Some parameter points are found to induce strong gravitational wave signals, which have the opportunity to be detected in future space-based interferometer experiments LISA, Taiji, and TianQin.


I. INTRODUCTION
It is conventionally believed that dark matter (DM) originates from thermal production at the early Universe [1][2][3]. Thus, the DM relic abundance would be determined by the annihilation cross section at the freeze-out epoch. The relic abundance observation suggests that the natural strength of the DM couplings to standard model (SM) particles should be close to the weak interaction strength. This motivates the worldwide establishment of various direct detection experiments searching for nuclear recoil signals induced by DM scattering. Nonetheless, no DM signal is robustly found in these experiments so far, leading to stringent constraints on the DM-nucleon scattering cross section [4][5][6]. Therefore, the thermal production paradigm faces a serious challenge.
Such a situation can be circumvented if one can effectively suppress DM-nucleon scattering at zero momentum transfer without reducing DM annihilation at the freeze-out epoch. An appealing approach to achieve this is provided by Higgs-portal pseudo-Nambu-Goldstone boson (pNGB) DM models [7][8][9][10][11][12][13][14][15][16][17][18][19][20], where the DM candidate is a pNGB protected by a global symmetry which is softly broken by quadratic mass terms. The pNGB nature makes the tree-level DM-nucleon scattering amplitude vanish in the zero momentum transfer limit [7]. Although loop corrections give rise to a nonzero scattering cross section, one-loop calculations show that near future direct detection experiments would not be able to probe the pNGB DM [8,9,19]. Therefore, other experimental approaches are crucial for exploring these models.
There are various ways to experimentally test DM models, including direct and indirect DM detection, collider searches, etc. The discovery of gravitational waves (GWs) by LIGO [21] provides a new path. DM fields could be relevant to strong first-order electroweak phase transitions (EWPTs) that produce detectable stochastic GW signals in the proposed future GW experiments . Such stochastic GWs typically peak around the mHz frequency band [50], to which ground-based laser interferometers are not sensitive. Nonetheless, future space-based GW interferometer plans, e.g., LISA [51], TianQin [52][53][54], Taiji [55,56], DECIGO [57,58], and BBO [59,60], are able to probe sub-Hz bands and look for the stochastic GW signals.
The minimal setup of the pNGB DM involves a complex scalar singlet with a global U(1) symmetry and a quadratic term that softly breaks U(1) into Z 2 [7]. The singlet and the SM Higgs doublet together could induce two-step phase transitions. Nevertheless, a study [12] showed that such phase transitions can only be of second order and impossible to produce stochastic GWs. Further studies tried to introduce extra terms to break the U(1) symmetry, e.g., the soft cubic terms [36], or the most general breaking terms [45]. These efforts successfully achieved first-order phase transitions (FOPTs) and stochastic GWs, but the essential merit of the vanishing tree-level DM-nucleon scattering in the zero momentum transfer limit is sacrificed.
It would be rather interesting if we can find out a pNGB DM setup that allows both the vanishing DM-nucleon scattering and detectable stochastic GW signals. Inspired by the notable GW signals from the strong FOPTs obtained in the two-Higgs-doublet models [61][62][63], we study the possibility of extending the minimal pNGB DM setup with an additional Higgs doublet [15], which is expected to involve more phase transition patterns. The corresponding phase transitions and stochastic GWs will be studied in this paper in detail.
In the following Sec. II, we briefly introduce the model and the particle masses. Existed experimental bounds are described in Sec. III. The effective potential at finite temperature is constructed in Sec. IV. We analyze key properties of the EWPT in Sec. V, which are relevant to the GW spectra discussed in Sec. VI. Sec. VII give numerical analyses of GW signals based on random parameter scans. We summarize the paper in Sec. VIII.

II. THE MODEL
In this section, we briefly describe the model we are interested in. More details can be found in Ref. [15]. This model involves two SU(2) L Higgs doublets Φ 1 and Φ 2 , both carrying hypercharge 1/2, and a complex scalar S which is a SM gauge singlet. The Lagrangian respects a global U(1) symmetry S → e iα S explicitly violated into a Z 2 symmetry S → −S by soft breaking quadratic terms. The U(1) symmetry is further spontaneously broken after S develops a vacuum expectation value (VEV). Then the imaginary part of S becomes a pNGB, acting as a DM candidate.
As in the simplified versions of the two-Higgs-doublet models [64], we assume that the scalar potential respects a Z 2 symmetry Φ 1 → −Φ 1 or Φ 2 → −Φ 2 which is only softly broken by quadratic terms. Moreover, CP conservation is assumed in the scalar sector, leading to only real coefficients. The potential satisfying these two assumptions and the global U(1) symmetry reads The U(1) soft breaking terms are further introduced in the potential. Thus, the total scalar potential is We can always make the soft breaking parameter m 2 S real and positive through a phase redefinition of S. Consequently, the potential respects a dark CP symmetry S → S * . For m 2 S > 0, the VEV of S developed must be real, and the dark CP symmetry remain unbroken, ensuring that the imaginary part of S acts as a stable DM candidate [7].
If the charged component of Φ 1 or Φ 2 gains a nonzero VEV, the photon would become massive, and the theory is unacceptable. If the neutral component of Φ 1 or Φ 2 develops an imaginary VEV, CP would be spontaneously broken. Detailed discussions on vacuum configurations and parameter relations in general two-Higgs-doublet models can be found in Ref. [65]. Here we are particularly interested in the case that only the neutral real parts of Φ 1 , Φ 2 , and S develop nonzero VEVs v 1 , v 2 , and v s , respectively. Thus, at the zero temperature, these scalar fields can be expanded as The potential is minimized at (v 1 , v 2 , v s ), leading to three stationary point conditions, where The VEVs contribute to a 3 × 3 mass-squared matrix for the CP -even neutral scalars ρ 1 , ρ 2 , and s. The eigenvalues m 2 h 1 , m 2 h 2 , and m 2 h 3 of this matrix are the masses squared for the mass eigenstates h 1 , h 2 , and h 3 , respectively. One of h i must behave as a SM-like Higgs boson with a mass of ∼ 125 GeV, satisfying the experimental observations. After rotations with the angle β, the CP -odd neutral scalars η 1 and η 2 are transformed into the mass eigenstates G 0 and a, while the charged scalars φ + 1 and φ + 2 are transformed into the mass eigenstates G + and H + . G 0 and G ± are the Nambu-Goldstone bosons eaten by the Z and W ± gauge bosons. a and H ± are extra Higgs bosons, whose masses squared are given by For m 2 S = 0, the neutral boson χ is a massless Nambu-Goldstone boson due to the global U(1) symmetry. The soft breaking terms endow the pNGB χ with a mass of Besides, χ only appears in pairs in the interaction terms, guaranteeing its stability to become a DM candidate. The pNGB feature also eliminates the tree-level χ-nucleon scattering amplitude in the zero momentum transfer limit without any parameter tuning [7,15]. Thus, this model is hardly constrained by DM direct detection experiments. The masses of the W and Z gauge bosons are given by where v ≡ v 2 1 + v 2 2 , and g and g denote the SU(2) L and U(1) Y gauge couplings, respectively. Thus, we observe that v is equivalent to the Higgs VEV in the SM and can be expressed as v = ( √ 2G F ) −1/2 , where G F is the Fermi constant. For the two Higgs doublets, four types of Yukawa couplings without tree-level flavorchanging neutral currents (FCNCs) can be constructed [64,66,67]. In this paper, we only focus on the type-I and type-II Yukawa couplings, whose Lagrangians are respectively given by where , and i, j = 1, 2, 3. The Yukawa coupling matricesỹ ij d andỹ ij u can be diagonalized by unitary matrices, which transform the gauge eigenstates u i and d i into the mass eigenstates u i and d i . We remark that due to the similarity of the Yukawa couplings in the quark sectors and the smallness of the ones in the leptonic sectors, many of the following analyses for the type-I (type-II) case can be cast to the lepton specific (flipped) case.

III. EXPERIMENTAL BOUNDS
In our analyses, we carry out random scans in the parameter space. The following 12 parameters are adopted as the free parameters: Each parameter point in the scans should be tested by existed experimental bounds. Firstly, we require that m 2 h i (i = 1, 2, 3), m 2 a , and m 2 H + should be positive to guarantee physical scalar masses. Moreover, in order to ensure that the scalar potential is bounded from below, the following conditions from copositivity criteria [68,69] should be satisfied: Furthermore, we require one of h i acting as the SM-like Higgs boson with a mass within the 3σ range of the measured value m h = 125.18 ± 0.16 GeV [70]. The numerical tool Lilith 2 [71,72] is used to test whether the SM-like Higgs boson is consistent with LHC run 1 and run 2 Higgs measurements from ATLAS and CMS. Parameter points excluded by the data at 95% confidence level (C.L.) are abandoned.
Although FCNCs have been forbidden at tree level, they can arise from loop corrections. In particular, the loops involving the charged Higgs boson H ± significantly contribute to the FCNC B-meson decays, depending on m H ± and tan β. The analysis by the Gfitter Group [73] shows that the strongest constraint on the type-I (type-II) Yukawa couplings comes from the measurement of the FCNC decay B d → µ + µ − (B s → µ + µ − and B → X s γ). We further reject the parameter points that are excluded at 95% C.L. by these flavor physics constraints.
Then we impose the constraints from DM phenomenology. We utilize FeynRules 2 [74] and the MadGraph5 aMC@NLO [75] plugin MadDM 3 [76] to calculate the prediction of the DM relic abundance. The observed value of the relic abundance from the Planck experiment is given by Ω DM h 2 = 0.1200 ± 0.0012 [77], where Ω DM is the ratio of the DM energy density to the critical density of the Universe and h is the Hubble constant in unit of 100 km s −1 Mpc −1 . Only the parameter points predicting the observed relic abundance are preserved. MadDM 3 is also used to compute the DM annihilation cross section σ ann v d with an average velocity of 2 × 10 −5 , which is corresponding to DM annihilation processes at dwarf spheroidal galaxies. The 95% C.L. upper limits on σ ann v in the bb channel from the γ-ray observations of dwarf galaxies by the Fermi-LAT satellite experiment and the MAGIC Cherenkov telescopes [78] are employed to test the parameter points.
Below, we study the effective potential, cosmological phase transitions, and gravitational waves for the parameter points surviving from all the experimental bounds above.

IV. EFFECTIVE POTENTIAL
In order to investigate the cosmological phase transitions in the model, we need to construct the effective potential. We assume that only the CP -even neutral scalar fields ρ 1 , ρ 2 , and s can develop VEVs in the cosmological history. The effective potential is then expressed as a function of the classical background fieldsρ 1 ,ρ 2 , ands.
The tree-level effective potential in terms of the classical fields derived from Eqs. (1) and (2) is Here, m 2 11 , m 2 22 , and m 2 S should be expressed as in Eqs. (6), (7), and (8), respectively. At zero temperature, the one-loop effective potential V 1 receives the Coleman-Weinberg terms [79] in the MS renormalization scheme [80], where the sum runs over all the particles i coupling to the classical fields, andm 2 i are the corresponding particle masses squared in terms of the classical fields. For the SM fermions, we only take into account the top and bottom quark contributions, and neglect all the other much smaller Yukawa couplings. Hence, all the particles we include in the calculations are Although the photon γ would not contribute to Eq. (24), its longitudinal mode can contribute to the daisy potential V D , which will be discussed below. µ is the renormalization scale. For transverse gauge bosons, C i = 1/2, while for longitudinal gauge bosons, scalar bosons and fermions, C i = 3/2. n i count the degrees of freedom of the particles, given by where the minus signs for n t and n b characterize the feature of fermion loops. The subscripts L and T denote the longitudinal and transverse polarizations of the gauge bosons.
In terms of the classical background fieldsρ 1 ,ρ 2 , ands, the elements of the symmetric mass-squared matrixM 2 h for the CP -even neutral scalar bosons are derived as The mass-squared matrix for the CP -odd neutral scalar bosons is withλ 345 ≡ λ 3 + λ 4 − λ 5 , while the mass-squared matrix for the charged scalar bosons is The eigenvalues of these matrices give the masses squared of the scalar bosons, i.e., The masses squared of the DM candidate χ is obtained as The mass squared of the W ± boson is given bỹ The mass-squared matrix of the B and W 3 gauge fields is After diagonalization, the masses squared of the Z boson and the photon arẽ For the type-I Yukawa couplings, the masses squared of the top and bottom quarks arẽ where the couplings y t = √ 2m t /v and y b = √ 2m b /v are defined the same as in the SM. For the type-II Yukawa couplings, the masses squared becomẽ Notice that loop corrections generally shift the values of the VEVs as well as the renormalized mass-squared matrix of the CP -even neutral scalar bosons. To keep them intact, we introduce the following counterterms [81,82], The nine counterterm coefficients are determined by the following nine equations at The masses of the Nambu-Goldstone bosons G 0 and G ± vanish at (ρ 1 ,ρ 2 ,s) = (v 1 , v 2 , v s ) in the Landau gauge, inducing logarithmic IR divergence terms in Eqs. (47) and (48) This problem is due to the ill-defined renormalized Higgs boson masses at p 2 = 0 with massless Nambu-Goldstone modes, and one can fix it by setting the momenta of the Higgs bosons on shell [83,84]. Similar problems exist in the effective potential with higher loops, and more details can be found in Refs. [85,86]. An approximate treatment is to give an IR cutoff Λ IR to the Nambu-Goldstone boson masses [81], i.e., to setm 2 Here, we take Λ IR to be the mass of the SM-like Higgs boson. Solving Eqs. (46)-(48), we obtain Thermal corrections to the effective potential are crucial for studying the EWPT. The one-loop finite-temperature effective potential [87] can be expressed as where T is the temperature and the functions J B and J F are defined as We also consider the daisy diagrams, which can be significant. The corressponding contribution to the effective potential can be estimated by [88,89] m 2 i are the field-dependent boson masses squared with thermal corrections in the hightemperature limit and can be derived bȳ whereM 2 X (ρ 1 ,ρ 2 ,s) represents the mass-squared matrices or masses squared in terms of the classical fields, and Π X (T ) denotes the thermal corrections toM 2 X . The subleading offdiagonal elements of Π X (T ) can be neglected [88,90]. The diagonal elements of Π X (T ) for the scalar bosons are derived as Here, y 1 and y 2 are the contributions from the Yukawa couplings. For the type-I and -II cases, they are given by Type II: The thermal corrections to the electroweak gauge bosons are Note that Π W 3 L and Π B L are the corrections to the diagonal elements ofM 2 W 3 ,B in Eq. (40).
Finally, we obtain the total effective potential 1 V. PHASE TRANSITIONS Based on the effective potential constructed in the previous section, we can study its evolution with temperature. At sufficiently high temperatures, the effective potential is minimized at the origin (ρ 1 ,ρ 2 ,s) = (0, 0, 0), implying the restoration of the electroweak gauge symmetry. As the Universe cools down, extra minima appear. In particular, if there are two coexisted minima separated by a high barrier, strong FOPT could take place and result in a stochastic GW background. We utilize the numerical package CosmoTransitions [92] to analyze the phase transitions. For each parameter point in the random scans, we verify whether or not the minimum (ρ 1 ,ρ 2 ,s) = (v 1 , v 2 , v s ) is the global one of the zero-temperature effective potential. The parameter points that fail this test are rejected. Then we use CosmoTransitions to trace the temperature evolution of the local minima.
In this model, the three classical CP -even neutral scalar fields would develop VEVs, typically leading to multi-step cosmological phase transitions. In Fig. 1, we demonstrate the temperature evolution of multiple phases for a benchmark point (BP), whose parameters can be found in the BP3 column of Table I  At T 460 GeV, the system stays at the red minimum with v 1 (T ), v 2 (T ), v s (T ) = (0, 0, 0), respecting the electroweak gauge symmetry. At T 460 GeV, a second-order phase transition occurs and the system turns into the green minimum, wheres develops a nonzero VEV. At T 148 GeV, the blue minimum appears, accompanied with a barrier that separates it from the green minimum. These two minima coexist till the zero temperature.
The effective potential at the blue minimum is higher than at the green minimum until the critical temperature T c 119 GeV. Below T c , the green minimum becomes a metastable state, i.e., a "false vacuum". The system finally undergoes a FOPT through quantum tunneling and turns into the blue minimum, or the "true vacuum". Such a FOPT nucleates bubbles, inside which the system is trapped at the true vacuum. In this FOPT, v 1 (T ) and v 2 (T ) increase from zero to O(100) GeV, while v s (T ) slightly decreases. At zero temperature, the true vacuum satisfies In our parameter scans, we usually find that v s (T ) does not evolve synchronously with v 1 (T ) and v 2 (T ), probably due to the less couplings of the singlet field to other fields compared with the two Higgs doublets. Typically, v s (T ) becomes nonzero much earlier than 1 Discussions on theoretical uncertainties in perturbative calculations of the effective potential can be found in Ref. [91]. the conventional EWPT epoch via a second-order or first-order phase transition.ρ 1 andρ 2 then gain VEVs in an subsequent phase transition, which could be a strong FOPT similar to those in the conventional two-Higgs-doublet models [93,94].
Below we discuss the dynamics of the FOPTs. The bubble nucleation rate per unit time and unit volume is given by [95,96]   action S 3 can be simplified to where r is the radius of the bubble. φ i (r) = ρ 1 (r),ρ 2 (r),s(r) is given by the bounce solution of the equations of motion with boundary conditions where φ false i is the field configuration of the false vacuum.
We present S 4 and S 3 /T as functions of the temperature for BP3 in Fig. 2(a), as well as the corresponding nucleation rate Γ in Fig. 2(b). We find that S 3 /T is the smaller one until temperatures below ∼ 5 GeV. The minimal of S 3 /T is reached at T ∼ 40 GeV. As the Universe cools down, below the critical temperature T c , the nucleation rate increases before the peak around T ∼ 40 GeV, and then decreases.
The bubbles are actually nucleated at the nucleation temperature T n , where the nucleation probability for a single bubble within a Hubble volume reaches O(1). Thus, T n can be estimated by [97] tn tc dt Γ where H is the Hubble rate, and t c and t n denote the critical and nucleation times, respectively. Note that the differential relation between the time t and the temperature T is dt = −(HT ) −1 dT in the radiation-dominated epoch.
Below the nucleation temperature T n , an increasing number of bubbles thrive and collide with each other. The maximum of bubble collisions that remarkably produces stochastic GWs is expected to be reached when percolation occurs [98]. In order to evaluate the percolation time t p , we need to estimate the fraction of space that still remains in the false vacuum at time t, which can be computed by [99,100] where a(t ) is the scale factor. r(t, t ) is the comoving radius of a bubble growing from t to t, given by where v w is the velocity of the bubble wall. For randomly distributed spherical bubbles with equal size in the three-dimensional space, the percolation threshold is reached when the fraction of space converted to the true vacuum, 1−P (t), increases to ∼ 0.29 [101,102]. Thus, the percolation time t p can be derived by requiring P (t p ) 0.71, with the corresponding temperature T p characterizing GW production from FOPTs [98,[103][104][105].
FOPTs are able to release latent heat from the vacuum energy, which drives the expansion of the bubbles and also converts into the thermal and bulk kinetic energies of the plasma [106][107][108]. The density of the released vacuum energy is given by [109] where φ true i is the field configuration of the true vacuum. It is useful to define a dimensionless strength parameter with ρ rad = π 2 g * T 4 /30 the radiation energy density in the plasma. g * is the effective relativistic degrees of freedom in the plasma.
The expansion of the bubbles depends on the interactions between the bubble walls and the plasma, analogous to chemical combustion in a relativistic fluid [106]. Hydrodynamic analyses show that bubble propagation have diverse modes, including Jouguet detonations, weak detonations, subsonic deflagrations, supersonic deflagrations (hybrid), and runway bubble walls [108]. Thus, it is difficult to completely work out the bubble wall velocity v w . For Jouguet detonations, the Chapman-Jouguet condition leads to a wall velocity of [106] v This is a typical assumption when evaluating GW signals.
Expanding the action S around the time t = t n or t = t p , we have where can be roughly understood as the inverse time duration of the phase transition [110]. For the electroweak FOPTs in which we are interested, the derivative dS/dT is positive at T = T , leading to positive β. In addition, S 3 /T is typically smaller than S 4 at T = T . In order to conveniently compare the phase transition time scale β −1 and the cosmological expansion time scale H −1 , we define a dimensionless quantitỹ Based on Eqs. (74) and (75), further calculations show that the nucleation and percolation temperatures T n and T p can be approximately determined by [111] For BP3, the nucleation temperature is T n 64 GeV, while the percolation temperature assuming v w = v CJ is slightly lower, T p 60 GeV, as denoted in Figs. 1 and 2.

VI. GRAVITATIONAL WAVE SPECTRA
Electroweak FOPTs could induce significant perturbations of the Friedmann-Robertson-Walker metric and produce stochastic GWs around the mHz band. Two key parameters relevant to the relic GW spectrum are α andβ evaluated at the time t * when GWs are produced. There are three coexisting GW sources at a FOPT, namely bubble collisions, sound waves, and magnetohydrodynamic (MHD) turbulence [112][113][114][115]. Denoting Ω GW to be the present GW energy density per logarithmic frequency interval divided by the critical density, we separate the contributions from the three sources as (85)

(a) Bubble collisions
The nucleated bubbles expand and finally collide with each other. Their collisions break the spherical symmetry and generate gravitational waves [110]. This process can be well described by the envelope approximation [110,116,117]. Numerical simulations for bubble collisions in the thermal plasma [107,118] show that the resulting GW spectrum at present can be approximated by where g * is evaluated at T = T * , the temperature corresponding to t = t * . The peak frequency of the spectrum can be modeled as [118] The redshift of the frequency has been taken into account by the factor h * = a(t * )H(t * ) a(t 0 ) = 1.65 × 10 −5 Hz T * 100 GeV g * 100 which is the inverse Hubble time at t = t * redshifted to today (t = t 0 ). The efficiency factor κ φ characterizes the fraction of the available vacuum energy converted into the gradient energy of the scalar fields.

(b) Sound waves
The explosive bubble expansion in the plasma induces a sound shell around the bubble wall. After the bubble collisions, the sound shells propagate into the fluid as sound waves, which become a significant GW source [119][120][121]. This source lasts until the sound waves are disrupted by the development of nonlinear shocks and turbulence [104,[121][122][123]. Therefore, the duration of the sound wave source can be determined by the nonlinearity timescale estimated as [123] where H * ≡ H(t * ) is the Hubble rate at t = t * and κ v is the fraction of the available vacuum energy converted into the kinetic energy of the fluid bulk motion. For Jouguet detonations, v w = v CJ , and κ v can be approximated by [108] In an expanding radiation-dominated Universe, the finite duration of the sound wave source leads to a suppression factor [124] Thus, the GW spectrum contributed by the sound waves is given by [121,124] where the peak frequency is estimated to be [121]

(c) MHD turbulence
Bubble collisions can stir up turbulence in the fluid, as the energy injection to the plasma results in an extremely high Reynolds number [107]. Since the plasma is fully ionized, the magnetic field, along with the velocity field, should be considered, leading to MHD turbulence [125]. It takes several Hubble times for the MHD turbulence to decay, and the stochastic GWs arise continuously during this period [126]. The corresponding GW spectrum can be fitted as [113,126] with Based on the suggestion from simulations, we optimistically set κ turb 0.1κ v [113,120].
In general, the contribution from the sound waves dominates in the GW spectrum [115]. Moreover, κ φ is typically negligible, except for runaway bubble walls [108,113,122]. Thus, we omit the contribution from the bubble collisions in the following calculations.  Type I  Type II  BP1  BP2  BP3  BP4 FIG. 3. Contours of the peak amplitudes of the GW spectra in the α-β −1 plane assuming Jouguet detonations. Purple and green points denote the parameter points for type-I and type-II Yukawa couplings, respectively. Four BPs are also indicated.
Note that positive λ 1 , λ 2 , and λ S are required to satisfy the bounded-from-below conditions (18). A negative v s would be totally equivalent to a positive one due to the Z 2 symmetry S → −S. Besides, a parameter point with tan β and m 2 12 is equivalent to one with − tan β and −m 2 12 , since the potential respects the Z 2 symmetry Φ 1 → −Φ 1 or Φ 2 → −Φ 2 expect for the soft breaking quadratic terms with m 2 12 . Thus, we can just take positive v s and tan β in the scans, while m 2 12 , λ 3 , λ 4 , λ 5 , κ 1 , and κ 2 can be either positive or negative. In addition, the v s range of 10 GeV to 1 TeV ensures that S has a VEV near the electroweak scale, and the interplay between S and the two Higgs doublets could be important.
The strength of the stochastic GW signals from the FOPT depend on α and β. Larger α implies a stronger FOPT, while smaller β corresponds to a longer FOPT time duration. Consequently, larger α andβ −1 lead to stronger GW signals, as implied in Eqs. (86), (92) and (94). For the surviving parameter points, we calculate the resulting values of α and β −1 , and then project the points in the α-β −1 plane, as presented in Fig. 3. The purple and green points are corresponding to type-I and type-II Yukawa couplings, respectively. The parameter points lie in the ranges of 10  FIG. 4. Peak amplitudes of the total GW spectra versus frequency for the parameter points with type-I (a) and type-II (b) Yukawa couplings assuming Jouguet detonations. Sensitivity curves for the future space-based GW interferometers LISA [51], Tianqin [54], Taiji [56], BBO [60], and ultimate DECIGO [58] are also plotted. The color axes denote the LISA signal-to-noise ratio SNR LISA for the parameter points with SNR LISA > 10. The gray points yield SNR LISA < 10.
We introduceΩ GW h 2 to denote the peak amplitudes of the GW spectra. The contours ofΩ GW h 2 are demonstrated in Fig. 3 and we can easily read off the GW signal strengths of the parameter points from the plot. The strongest GW signal we find reaches up tô Ω GW h 2 ∼ 10 −11 . Fig. 4 illustratesΩ GW h 2 versus the peak frequency f for the parameter points. For comparison, we also plot the sensitivity curves for the future space-based interferometers LISA [51], TianQin [54], Taiji [56], BBO [60], and DECIGO [58]. Some of the curves are converted from the sensitivity on amplitude spectral density or characteristic strain. The conversions of the related quantities can be found in, e.g., Ref. [127]. The DECIGO curve we adopt here is the ultimate sensitivity that is only limited by quantum noises, and it can be regarded as an observational limitation [58].
The GWs produced by FOPTs become an isotropic and stochastic background in the present Universe. The detectability of the GW signals in the space-based interferometers increases with the practical observation time T . The signal-to-noise ratio can be defined as [113,128] SNR ≡ T fmax f min where Ω sens (f ) is the sensitivity of the experiment. Below we take the practical observation time T = 9.46 × 10 7 s (3 years) for LISA [114], Taiji, and TianQin. The signal Ω GW (f ) can be detected if the corresponding SNR is larger than a signal-to-noise ratio threshold SNR thr . For the six (four) link configuration of LISA, the threshold is SNR thr = 10 (50) [113]. We find that some parameter points yield the LISA signal-to-noise ratio SNR LISA > 10 and could be probed by LISA. We denote SNR LISA for them as the color axes in Fig. 4, with the remaining gray points corresponding to SNR LISA < 10. The next-generation plans aiming at f ∼ O(0.1) Hz, like BBO and DECIGO, may probe much more parameter points.
For a closer look at the results, we choose four benchmark points, whose detailed information is listed in Table I. BP1 and BP2 (BP3 and BP4) correspond to the type-I (type-II) Yukawa couplings. In these BPs, the masses of the Higgs bosons h 1,2,3 , a, and H ± are all below 1 TeV, while the mass of the DM candidate is less than 140 GeV. The SM-like Higgs boson is h 2 in BP2, while it is h 1 in the rest BPs. The DM annihilation cross sections σ ann v d at dwarf galaxies predicted by the BPs are below 2 × 10 −26 cm 3 /s, beyond the reach of Fermi-LAT and MAGIC [78].
For the four BPs, percolation of the FOPT occurs in 47 GeV T p 75 GeV, with α ranging from 0.16 to 0.35 andβ −1 ranging from 4 × 10 −3 to 2.2 × 10 −2 . Assuming Jouguet detonations, we derive the GW spectra for these BPs, as presented in Fig. 5. The BPs are also indicated in Figs. 3 and 4. We find that the GW signal strengths decrease according to the order of BP4, BP1, BP3, and BP2, reflecting the descending order of α. In Table I, we also list the signal-to-noise ratios SNR LISA , SNR Taiji , and SNR TianQin for the LISA, Taiji, and TianQin experiments, respectively. LISA and Taiji look promising to detect all BPs, while TianQin may probe BP4 with a sightly longer observation time.
In order to show the most important flavor constraints mentioned in Sec. III, we plot our parameter points confronting the FCNC bounds. We have adopted the data from the Gfitter global fit [73] to reject parameter points. In Fig. 6, the parameter points for the two types of Yukawa couplings are projected in the tan β-m ± H plane, with the color axis indicating the peak amplitude of the GW spectrum,Ω GW h 2 .
For the type-I case in Fig. 6(a), the most stringent FCNC bound comes from the LHCb and CMS measurements of B d → µ + µ − [129,130], excluding a region with tan β 3. For the type-II case in Fig. 6(b), the bounds from the observations of B → X s γ [131][132][133] and B s → µ + µ − [129,130] exclude a region with m H ± 750 GeV. Thus, the FCNC constraints remove small tan β and light H ± for type-I and type-II Yukawa couplings, respectively. We remark that strong GW signals typically favor small m H ± . There are many parameter points with m H ± 600 GeV in the type-I case leading to largeΩ GW h 2 . In the type-II case, the m H ± 600 GeV region is basically excluded by the FCNC constraints, and hence it is more    difficult to achieve strong GW signals. In Fig. 7, we project the parameter points in the m χ -σ ann v d plane. Although all these parameter points are required to predict the observed relic abundance, σ ann v d can deviate from the canonical annihilation cross section 3×10 −26 cm 3 /s, due to the velocity dependence of σ ann v . One reason is that some annihilation channels kinematically forbidden at low velocities could be opened at the freeze-out epoch, and another reason is the resonance effects [134]. The pile-up of points around m χ ∼ 78 GeV is mainly related to the annihilation channel χχ → W + W − .
In the above numerical analyses, we have assumed the bubble propagation mode to be Jouguet detonations with bubble wall velocity v w = v CJ . Below, we study the effect of various bubble propagation modes. In general, the dependence of κ v on v w and α can be found in Ref. [108]. The sound speed c s in the relativistic plasma is very close to 1/ √ 3. For v w c s , v w = c s , and v w = 1, the vacuum energy fraction converted into the bulk kinetic energy of the fluid κ v has the following analytic approximations, based on fit.
According to these expressions, we derive the GW spectra for BP4 assuming the bubble wall velocity v w = 0.05, 0.2, 0.72, and 1. These spectra are demonstrated in Fig. 8, along with the previously obtained BP4 spectrum for v w = v CJ = 0.87. Compared to Jouguet detonations with v w = v CJ , supersonic deflagrations with v w = 0.72 lead to a stronger GW signal, which could be properly tested by TianQin with SNR TianQin = 15.8. On the other hand, subsonic deflagrations with v w = 0.05 and 0.2 give much weaker GW signals.

VIII. SUMMARY
In this paper, we have studied the stochastic GW signals from electroweak FOPTs in the model comprising the pNGB dark matter framework and two Higgs doublets with type-I or type-II Yukawa couplings. The DM candidate is a pNGB whose tree-level scattering off nucleons vanishes at zero momentum transfer, evading the constraints from direct detection experiments. The three scalar fields in the model have nonzero VEVs at zero temperature, which should be developed from EWPTs at the early Universe. If such EWPTs are of strongly first order, stochastic GWs could be effectively produced. The effective potential has been carefully constructed with one-loop corrections at zero temperature as well as thermal corrections, allowing us to carry out accurate analyses on EWPTs.
We have performed random scans in the 12-dimensional parameter space, taking into account the constraints from bounded from below conditions, LHC run 1 and run 2 measurements of the 125 GeV Higgs boson, FCNC B-meson decays, the Planck observation of the DM relic abundance, and the γ-ray observations of dwarf galaxies by Fermi-LAT and MAGIC. The surviving parameter points are also required to induce an electroweak FOPT. We have further analyzed the characteristic temperatures, the phase transition strength, and the characteristic time duration of the FOPT. Based on such information, the resulting relic GW spectra from sound waves and MHD turbulence have been evaluated.
Assuming that the bubble propagation mode is Jouguet detonations with v w = v CJ , we have found that the FOPTs of some parameter points could induce peak amplitudes of the GW spectra around 10 −13 -10 −11 , which could be well detected by the future space-based GW interferometers LISA and Taiji. The next-generation GW interferometers BBO and DECIGO are capable of probing much more parameter points. We have noticed that a lighter charged Higgs boson H ± in this model is more probable to induce a strong GW signal. Since the FCNC constraints on m H ± for type-II Yukawa couplings at large tan β are more stringent than those for type-I Yukawa couplings, the type-I case typically leads to stronger GW signals.
We have also investigated the effects of different bubble propagation modes with several values of the bubble wall velocity v w . For the benchmark point BP4, supersonic deflagrations with v w = 0.72 can induce a stronger GW signal than Jouguet detonations with v w = v CJ . In this optimistic case, BP4 could be well tested by LISA, Taiji, and TianQin. Detonations with v w = 1 lead to a slightly weaker GW signal, while subsonic deflagrations with v w = 0.2 and 0.05 result in much weaker signals.