Electroweak Phase Transition with an SU(2) Dark Sector

We consider a non-Abelian dark SU(2)$_{\rm D}$ model where the dark sector couples to the Standard Model (SM) through a Higgs portal. We investigate two different scenarios of the dark sector scalars with $Z_2$ symmetry, with Higgs portal interactions that can introduce mixing between the SM Higgs boson and the SM singlet scalars in the dark sector. We utilize the existing collider results of the Higgs signal rate, direct heavy Higgs searches, and electroweak precision observables to constrain the model parameters. The $\text{SU(2)}_{\text{D}}$ partially breaks into $\text{U(1)}_{\text{D}}$ gauge group by the scalar sector. The resulting two stable massive dark gauge bosons and pseudo-Goldstone bosons can be viable cold dark matter candidates, while the massless gauge boson from the unbroken $\text{U(1)}_{\text{D}}$ subgroup is a dark radiation and can introduce long-range attractive dark matter (DM) self-interaction, which can alleviate the small-scale structure issues. We study in detail the pattern of strong first-order phase transition and gravitational wave (GW) production triggered by the dark sector symmetry breaking, and further evaluate the signal-to-noise ratio for several proposed space interferometer missions. We conclude that the rich physics in the dark sector may be observable with the current and future measurements at colliders, DM experiments, and GW interferometers.


Introduction
The milestone discovery of the Higgs boson predicted in the Standard Model (SM) at the CERN Large Hadron Collider (LHC) has deepened our understanding of nature at the shortest distances, and in the same time sharpened our questions about the Universe. One of the most pressing mysteries in contemporary particle physics and cosmology is the nature of the dark matter (DM). There is mounting evidence for the existence of DM through its gravitational effects. However, the null results of the last fifty years of searches challenge the most theoretically attractive candidates, namely, the standard weakly interacting massive particles (WIMPs), that are charged under the SM weak interactions (see Ref. [1] for review). On the other hand, it is quite conceivable that the DM particles live in a dark sector that are not charged under the SM gauge group. Furthermore, the dark sector may have a rich particle spectrum, leading to other observable consequences [2]. A massless dark gauge field, dubbed as the dark radiation (DR), is one of the quite interesting extensions that could help to alleviate the tension between Planck and HST measurements of the Hubble constant [3]. DM-DR interactions and DM self-interactions can provide solutions to the small-scale structure problems which challenge the cold dark matter (CDM) paradigm [4][5][6].
In this paper, we would like to explore the potentially observable effects beyond the gravitational interactions from a hypothetical dark sector. We assume that the dark sector interacts with the SM particles only through the Higgs portal [7]. An immediate consequence of this would be the modification of the Higgs boson properties that will be probed in the on-going and future high energy experiments [8,9]. The DM searches from the direct and indirect detection experiments will provide additional tests for the theory [1]. Perhaps, an even more significant impact would be on the nature of the electroweak phase transition (EWPT) at the early Universe (see, e.g., [10][11][12] for recent reviews), which could shed light on another profound mystery: the origin of baryon asymmetry in the Universe. Indeed, one of the best-motivated solutions to this mystery is the electroweak baryogenesis (EWBG) [13][14][15][16] (see also [17,18] for pedagogical introductions). For a successful generation of the baryon asymmetry during the EWPT, all of the three Sakharov conditions [19] have to be satisfied. One of the three Sakharov conditions is to assure a strong first-order phase transition (FOPT), that is absent within the minimal SM, but could be achieved by the Higgs portal to a sector beyond the SM. It is important to note that many well-motivated extensions of the SM predict gravitational wave (GW) signals through a strong FOPT, that are potentially detectable at LIGO and future LISA-like space-based GW detectors.
Given the rich physics associated with a dark sector, there have been significant activities in the literature dealing with many different aspects of the theory and phenomenology. We collect a broad class of sample models in Table 1. We classify them according to the particle charges under the SM gauge interactions, the dark gauge groups (D), and their state representations in the scalar sector, and the gauge symmetry breaking patterns. Those extensions with new states charged under the SM gauge group will yield substantial observable effects if kinematically accessible. Prominent examples include the extended Higgs sectors and supersymmetric theories. Those extensions uncharged under the SM gauge group will be characterized as dark sector. In the dark sector, both Abelian (U(1) D ) and non-Abelian (SU(2) D , SU(3) D ) gauge sectors have been studied with different symmetry breaking patterns induced by various scalar scenarios, as listed in Table 1. In particular, we examine their physical implications of FOPT, detectable Models Strong 1 st order GW signal Cold DM Dark Radiation and phase transition small scale structure
Building upon the existing literature, in this paper, we will focus on a dark SU(2) D model uncharged under the SM gauge group. Some early exploration and the phenomenology associated with the model have been examined [60][61][62][63][64][65][66][67][68][69][70]. The previous works mainly focused on the DM studies. In this work, we will study the EWPT and GW with this well-motivated DM model. In this class of models, it remains largely unconstrained on the choice of the dark scalar sector. With just one real scalar triplet, we could achieve a FOPT at the early Universe by transitioning from an electroweak symmetric vacuum that breaks the SU(2) D symmetry to an electroweak broken vacuum that preserves the SU(2) D symmetry [76]. As such, all the dark sector particles would remain massless, and there would be no cold DM candidate in this simplest scenario. Alternatively, we would like to explore the following two cases to facilitate a strong FOPT in the early Universe and to have viable cold DM candidates 1. one real scalar triplet and one real scalar singlet; 2. two real scalar triplets.
For both cases, at zero temperature, only one scalar triplet gets a nonzero vacuum expectation value (VEV) and partially breaks the SU(2) D into U(1) D . The massless vector gauge boson associated with the unbroken U(1) D symmetry can serve as a dark radiation (DR). The other two massive gauge bosons associated with the symmetry breaking are our vector DM candidates. Due to the presence of the non-Abelian gauge boson couplings, the DM-DR and DM-DM interactions can be naturally introduced. The other scalar triplet or singlet can develop a non-zero VEV at a finite temperature and can thus trigger a strong FOPT, besides providing the scalar DM candidates. We also list the main features of our model in the last row of Table. 1.
The rest of the paper is organized as follows. In section 2, we introduce our model and particle spectrum, with the phenomenological constraints presented in section 3 and DM phenomenology in section 4. In section 5, we perform the study of EWPT and the GWs spectrum with two benchmark points (BMs) as shown in Table 2. We summarize in section 6.

Theoretical framework
In addition to the SM, we include a non-Abelian SU(2) D dark sector. We consider two scenarios for the dark scalar sector, a real singlet plus a real triplet (ST), or two real triplets (TT) under the dark gauge group SU(2) D : We assume that the dark sector does not carry SM charges but rather interacts with the SM particles through the Higgs portal interactions. Therefore, the Lagrangian of the model consists of three parts µW c ν is the dark gauge field strength tensor; D µ = ∂ µ − igT aW a µ is the covariant derivative in the dark sector with T a being the SU(2) D generators, which is given in the 3-dimensional representation by (2. 6) and , being the SM Higgs doublet. The most general renormalizable hidden sector potential with an assumed Z 2 symmetry is given by where λ 4 = 0 in the ST model. In principle, there can be cubic terms for the singlet scalar, which can change the phase transition dramatically. However, we will not consider breaking the Z 2 symmetry in this work. 1 In our phenomenological analyses in the following sections, we choose v 1 = 0 at the zero temperature. An important consequence of this choice is to leave the dark U(1) D unbroken so that there will be a massless dark gauge field, DR, which would have observational implications.

Mass spectrum
With the choice of v 1 = 0, the SM Higgs boson mixes only with the SU(2) D dark scalar ϕ 2 . In the TT scenario, the mass terms for the scalar bosons are where h = {h 0 , ϕ 2 } are two neutral scalars with the mass matrix (2.9) and m 2 ω ± = 1 2 (λ 3 v 2 2 + 2m 2 11 + λ H11 v 2 h ) is the mass of the SU(2) D charged scalars. The mass for another neutral scalar ω 2 is m 2 (2.10) In the ST scenario, there is only one massive scalar with mass Please note that the sign ± refers to the dark SU(2) D charge. The neutral scalars h 0 and ϕ 2 are mixed. The mass eigenstates h = {h 1 , h 2 } can be obtained from a rotation on h The rotation matrix can be parametrized by one mixing angle θ as The mass eigenvalues are (2.14) Here and henceforth, we identify h 1 as the SM-like Higgs boson with m h 1 = 125 GeV, and h 2 is a heavier scalar in the model. The scalar fields ϕ 1 and ϕ 3 are the Nambu-Goldstone (NG) bosons absorbed by two of the SU(2) D gauge bosonsW 1 andW 3 . The mass terms of dark gauge bosons are contained in ( andW 2 remains massless.

Interactions
The interactions between the SM and the dark sector are generated through the Higgs portal as in Eq. (2.4), specifically where the scalar couplings are given in terms of the mixing angle and the other model parameters In the ST scenario, c i , c ij , and λ 4 are zero. The above interactions govern the phenomenology relevant for the potential experimental observations, such as the Higgs properties, the DM relic density and direct detections, and EWPT at the early Universe, as we will explore in the following sections.

Phenomenological constraints
The scalar potential of the model is The two minima conditions ∂V S ∂h 0 = 0 and ∂V S ∂ϕ 2 = 0 evaluated at the VEVs are The mass parameters m H and m 22 can be solved by using these two minima conditions In the TT model as described in the last section, there are fourteen parameters g, v h , v 1 , v 2 , m 2 H , m 2 11 , m 2 22 , λ H , λ H11 , λ H22 , λ 1 , λ 2 , λ 3 , λ 4 .
By applying the two extrema conditions in Eqs. (3.2) and (3.3) for the scalar potential and v 1 = 0, we can get rid of three parameters. Adopting the SM values m h 1 = 125 GeV, v h = 246 GeV, we are left with nine independent parameters, which can be chosen as sin θ,g, mW + , m h 2 , m ω + , m ω 2 , λ 1 , λ H11 , λ 3 . In the ST model, we have one less free parameter as m ω + and m ω 2 are replaced by one parameter m ω . We wish to have observable imprints from the dark sector in the current and future experiments. We thus take the SU(2) D symmetry breaking not too far from the electroweak scale in the SM, and vary the mass of the second Higgs boson m h 2 in the range of 200 GeV−1 TeV. We will not consider m h 2 > 1 TeV, as the perturbative GW calculations are not reliable. We examine the possible bounds on the other model parameters from the existing experiments in the following sessions. For the purpose of illustration, we choose two benchmark points (BMs) for the input parameters as shown in Table 2

Vacuum stability
A stable physical vacuum has to be bounded from below keeping the scalar fields from running away. The behavior of the scalar potential is dominant by the quartic part when the field strength approaches infinity. The conditions of vacuum stability are given in Ref. [78,79]. Following their procedure, we find the following conditions

Partial wave unitarity
The scattering amplitudes for spin-less 2 → 2 processes can be decomposed into a sum over the partial waves a j as where P j (cos α) are the Legendre polynomials in terms of the scattering angle α. The perturbative unitarity requires Im(a j ) = |a j | 2 , which implies We will adopt the second condition as it turns out to be more constraint. The s-wave amplitude can be computed by For a spin-less 2 → 2 elastic scattering process, the unitarity bound can be rephrased as Owing to the Goldstone-boson equivalence theorem, the scattering of the longitudinal gauge bosons can be approximated by the pseudo-Goldstone boson scattering in the high-energy limit. Given the fact that the high energy scattering is dominated by the four-scalar contact interactions, we only need to evaluate the quartic or bi-quadratic terms. There are ten scalar fields in the TT scenario, namely, ω i (i = 1 to 3), ϕ j (j = 1 to 3), G k (k = 0 to 2), and h 0 . So there are 55 pair combinations and 1540 scattering channels. An additional symmetric factor 1/ √ 2 needs to be included for each pair of identical particles in the initial or final states. The unitarity bounds from scattering amplitude matrix A 55×55 are 14) Similarly, for the ST case, there are a total of eight scalar fields and therefore 36 pair combinations. The unitarity bounds from scattering amplitude matrix A 36×36 are The horizontal purple line is from the Higgs signal rate measurement [83]. The yellow shaded region shows the upper bound from the direct searches for the heavy Higgs at LEP and LHC ( √ s = 7 TeV) [84]. The blue (red) shaded regions are excluded by the LHC di-boson searches with VBF (ggF) channels. The blue and red dashed lines correspond to the HL-LHC projection for these two channels, respectively [85]. The grey shaded area labelled by W mass, and the area above the brown dashed line labelled by S, T, U are excluded by the electroweak precision observables [80].

Electroweak precision observables
Quantum corrections to the W boson mass [80] and the electroweak oblique parameters [81], from the mixing between SM Higgs and the dark massive eigenstates, can put constraints on the model parameters sin θ and m h 2 . The bound from W boson mass constraint, which is shown by the gray shaded region in Fig. 1, turns out to be more stringent than that from the oblique parameters [80,82]. The bound from oblique parameters are shown by the dashed brown line in Fig. 1 for comparison.

Higgs phenomenology
The scalar state h 0 mixes with ϕ 2 after the electroweak symmetry breaking. We identify that the lighter mass eigenstate h 1 is the observed SM-like Higgs boson with a mass of 125 GeV. The couplings of the physical scalars h 1 and h 2 to the SM particles are The SM-like Higgs boson coupling to the SM particles are modified by a universal factor cos θ.
The relevant Higgs self-interactions in the scalar sector are where v 2 = mW + /g. These couplings are important for the DM annihilation at the early Universe through the Higgs portal. The Higgs phenomenology at colliders is similar to that of one real singlet scalar extension of the SM, which has been extensively studied (see [41,45,86,87] and references therein). The most relevant parameters are the mixing angle θ and the mass of the second Higgs m h 2 as shown in Eq. (3.19). The current bounds on sin θ and m h 2 from the Higgs phenomenology are shown in Fig. 1. We will discuss the details of each bound in the following subsections.

Higgs invisible decay
In the case that DM masses are larger than the half of the Higgs boson mass, the invisible decay of the Higgs boson is to the DRW 2 through the SU(2) D charged scalar and gauge bosons loops as shown in Fig. 2. The decay width through dark gauge bosons can be calculated as In the limit m h 1 m ω + , the decay width through dark scalars can be calculated as where c 1 is the coupling of vertex h 1 ω + ω − given in Eq. (2.18). The Higgs invisible decay width for our benchmark points shown in Table 2 are (3.25) which are dominated by the last two diagrams in Fig. 2. The Higgs invisible decay is highly suppressed by the small mixing angle, dark-sector gauge coupling, and the one-loop suppression. The branching fractions of the invisible Higgs decay are far beyond the reach of current and future experiments. Figure 2: Feynman diagrams for the Higgs invisible decay to the dark radiation.

Higgs coupling measurements
Higgs couplings with SM particles have been measured with good precisions at the LHC. The Higgs signal strength is defined as [88] Therefore, the signal strength can be written as As we learned from the previous section, Γ DS h 1 are highly suppressed, as the SM-like Higgs h 1 can only decay to DR through one-loop diagrams in Fig. 2. The signal strength simply scales as cos 2 θ. The bound on the mixing angle, from the Higgs couplings measurement by ATLAS [83], is | sin θ| 0.35, which is shown by the purple line in Fig. 1.
Of special interest is the SM-like Higgs triple coupling κ 111 as in Eq. (3.20) because of its sensitivity to the BSM new physics and its crucial role in EWPT. We write the derivation from the SM prediction as We depict the resultant deviation of ∆κ 3 in the v 2 -sin θ plane in Fig. 3 by the gray solid lines. For most of the viable parameter space, the magnitude of ∆κ 3 is less than 25%. We also mark the predictions of our benchmark points BM1 for about −10% by the red-cross and MB2 for about −2% by the blue-star, respectively. The achievable sensitivity to probe ∆κ 3 in the future collider experiments has been extensively studied. While the HL-LHC will only have a moderate sensitivity to κ 3 [89,90], future improvements are highly anticipated, reaching a 1σ sensitivity of 13% at a 1-TeV ILC [91] and 10% at CLIC [92], and 2σ sensitivity of 5% at FCC hh /SPPC [93], 2% at a multi-TeV muon collider [94]. The precision measurement for κ 3 would provide important indirect test of the model as well as BSM theories in general.

Direct searches for the heavy Higgs boson
The heavy Higgs boson in the model, at the high-energy colliders can put strong constraints in this scenario. Heavy Higgs h 2 mainly decay to heavy particles when they are kinematically allowed, such as bb, top quarks, massive gauge bosons, and the dark gauge bosons. The branching fractions of the heavy Higgs decay versus m h 2 are shown in Fig. 4, where the other parameters are fixed as BM1 in Table 2 for illustration. The heavy Higgs decay channels are to di-bosons W W + ZZ until the threshold for W +W − is open, as shown in Fig. 4. The LHC di-boson resonance search in gluon-gluon fusion (ggF) and vector boson fusion (VBF) [85] can put strong bounds on the mixing angle θ and the heavy Higgs mass m h 2 . We evaluate the resonance production rate as The bounds on the plane in m h 2 -sin θ with v 2 = 1000 GeV are shown by the red (ggF) and blue (VBF) shaded regions in Fig. 1. The dashed lines with the same color scheme are the projected limit from HL-LHC with 3 ab −1 integrated luminosity, obtained by rescaling the current bounds by the square root of luminosity ratio 3000/36.1. For the mass below 350 GeV, we adopted the bounds provided in Ref. [84] from a combination of various decay channels at LEP and LHC with √ s = 7 TeV. The bounds are shown by the yellow shaded region.

Dark radiation and dark matter phenomenology
The dark sector in our model possesses rich phenomenology. There are two self-interacting vector dark matter candidatesW ± . The massless stateW 2 is the DR. In addition, there is one scalar DM ω in the ST scenario, or there are three self-interacting scalar DM candidates ω ± and ω 2 in the TT scenario. The DM interactions with SM particles are through the mixing between  Table 2. sector among scalar DM, vector DM, and DR, which are from the dark gauge couplings and the scalar potential. For simplicity we decoupled the scalar DM by assuming the mass hierarchy to be mW + m ω in ST, or mW + m ω + m ω 2 in TT.

Dark radiation
The massless DRW 2 associated with the unbroken U(1) D can contribute to the energy density of the Universe, regulating the Universe expansion rate. In the radiation-dominated era, the expansion rate of the Universe depends on the relativistic energy density where g * is the total relativistic degrees of freedom defined as where the coefficients are C i =1 (7/8) for bosons (fermions), and g i is the internal degrees of freedom for particle i.W 2 can contributes to g * and it is conventional to define this extra energy density by whereT is the dark sector's temperature, T ν is the SM neutrinos' temperature. After neutrinos decouple from the thermal bath, the ratioT /T ν is fixed as they evolve in the same way. We thus evaluate this temperature ratio at the epoch of neutrino decoupling. Before the DM decouples, the dark sector and visible sector are in thermal equilibrium,T dec,χ = T dec,χ . After decoupling of DM, the dark sector and visible sector lost thermal contact, the entropy is conserved in each sector separately. So we have [95] g DS * s (T dec,ν )T 3 where g * s is the relativistic degrees of freedom for entropy Currently, the strongest bounds on N eff come from the Planck satellite [96,97] which measured N eff = 2.99 ± 0.17 including baryon acoustic oscillation data. The projected limit of CMB Stage IV experiments is ∆N eff = 0.03 [98], which has sufficient sensitivity to explore this scenario.

Relic density
The observed value of the DM relic density Ω obs h 2 0.12 inferred by the Planck collaboration from the analysis of Cosmic Microwave Background (CMB) [99]. The vector DM candidates W ± and scalar DM candidates ω (ω ± and ω 2 ) in the ST (TT) model can account for the DM relic density we observed today. 2 By solving the Boltzmann equation in the standard freeze-out scenario, the relic density of our DM candidates can be estimated by [100] Ω DM h 2 = 1.07 × 10 9 x f GeV −1 (g * S / √ g * )M pl σv rel , (4.7) where x f ≡ m χ /T f , which can be estimated by Here g * (g * S ) is the effective degree of freedom in energy density (entropy) at freeze-out defined in Eq. (4.2) ((4.5)). We evaluate the s-wave annihilation cross section at the leading order [101] σv rel = 1 32π (4.9) The attractive longe-range force between the vector DMW ± introduced by the exchange of massless DRW 2 can increase the annihilation cross section, which is the so-called Sommerfeld enhancement. The Sommerfeld factor is given by [102] S =α π v (4.10) When the DM freezes out, we can safely ignore the effects of the Sommerfeld enhancement in this work for the relic density calculation. We calculated the annihilation cross section of the process The representative Feynman diagrams for the vector DMW ± pair annihilation are shown in Fig. 5. Scalar DM pair annihilations have similar diagrams. Since we choose m ω ± , m ω 2 mW ± , scalar DM candidates ω ± and ω 2 will be decoupled much earlier than vector DMW ± . The scalar DM states in the TT model annihilate dominantly into the vector DMW ± . While, in the ST model, the scalar DM annihilation channel is dominated by ωω → h 2 h 2 as it does not carry any charge. At the decoupling of ω ± and ω 2 , nW ± = n eq W ± . Therefore, including the DM self-interacting processes can further reduce the relic density of ω ± and ω 2 . The number densities of ω ± and ω 2 are much less thanW ± at the decoupling ofW ± . Therefore, we ignore the processes ω + ω − →W +W − and ω 2 ω 2 →W +W − when we evaluate the number density ofW ± . The vector DM mainly annihilates into the DRW 2 except in the resonance region mW ± ≈ m h 2 /2. The branching fractions to a specific final state from an initial state annihilation of both vector and scalar DM pairs are shown in Fig. 6.
The relic densities for some benchmark points are shown in the left panel of Fig. 7 as functions of heavy Higgs mass m h 2 . The vector DM relic density is highly suppressed at the resonance  Table 2. region. The scalar DM contributions to the total relic density are negligible. The dashed green lines are the scalar DM from the ST scenario, which mostly overlaps with ω ± as they have the same masses and similar annihilation channel as shown in Fig. 6. We require the DM not to be overly produced Ω DM h 2 0.12. The dashed horizontal line in the left panel of Fig. 7 indicates the current relic density bound from PLANCK. In the resonance region mW ± ≈ m h 2 /2, the annihilation cross sections via an s-channel h 2 are enhanced, and the relic density is much less than the observed value. Away from the resonant region,W ± could be adequate as a CDM candidate.

Direct detection
The null results of direct detection experiments can set strong bounds on our dark sector parameter space. In this model, the DM candidates χ couple to the dark scalar ϕ 2 . ϕ 2 couples to the SM particles through the Higgs portal. The dominant contributions to the spin-independent (SI) scattering cross section come from the exchange of the SM-like Higgs bosons h 1 and the heavy Higgs bosons h 2 . The effective interactions of DM (χ =W ± , ω ± , ω 2 , ω) with light quarks and gluons are given as [1] L eff q,g = q=u,d,s  Table 2. where G a µν is the field strength tensor of gluon and α s is the strong coupling constant. f χ q is the effective couplings between DM χ and light quarks, which, in our model, are ). (4.14) The coupling between DM and gluon comes from the effective coupling after integrating-out of heavy quarks The interactions between DM and nucleon can be evaluated by using the nucleon matrix elements where f N T q and f N T G are the mass-fraction parameters of the quarks and the gluon in the nucleon N ,respectively. In our numerical calculations, we adopt f p T d = 0.0191, f p T u = 0.0153, f p T s = 0.0447, and f p T G ≡ 1 − q=u,d,s f p T q = 0.925 [103]. The effective interactions of DM and nucleon can be expressed as where the effective coupling f N can be calculated by The SI cross section of DM with nucleon can be calculated with [104] 19) where m N is the mass of nucleon and m χ is the mass of DM candidate. To derive the experimental upper bound, we scale the SI cross sections with the density fractions The XENON1T [105] and the SI cross sections are shown in the right panel of Fig. 7. In the resonance region mW ± ≈ m h 2 /2, the relic density is much less than the observed value, hence the direct detection bound can be easily evaded. Away from the resonant region however,W ± could lead to a detectable cross section.

Dark matter self-interactions
The collision-less and cold DM can successfully describe the large scale structure of the Universe [106]. There are, however, some challenges for the cold and collision-less DM model at the small-scale (see Ref. [107] for a review). Rather than going to the warm DM scenario, there are generally two mechanisms which can alleviate the CDM challenges: (i) DM-DR interactions [5]; (ii) DM self-interactions [4]. In our model, the leading DM self-interaction is mediated by the massless DR. This scenario has been studied carefully in Ref. [108,109]. The most relevant DM self-interactions are through t/u-channel mediated by the massless DR. The differential cross section of t-and u-channel in the center-of-mass (CM) frame is leading to σ ∼ πα 2 /(m 2W ± v 4 r ), where v r is the relative velocity of the two colliding DM particles in the CM frame. The cross sections of the DM self-interactions quickly drop at higher velocities to evade impacts on the large scale structure, hence, maintain the effective collision-less descriptions. From the observed ellipticity of galactic DM halos [108,109], a bound on the dark gauge couplings can be estimated as This constraint can potentially be overly strong and depends on the assumptions of DM relic density [109]. The constraints from the Bullet Cluster are much weaker [108,109]. To solve the small-scale structure problems, we need σ/mW ∼ 0.1 − 10 cm 2 /g at dwarf galaxies [64,110], which gives DM can also interact with themselves through four-gauge-boson contact and s-channel interactions. The cross sections of contact interactions are σ ∼ πα 2 /m 2W ± ; the s-channel cross sections are σ ∼ πα 2 v 4 r /m 2W ± . Therefore they are irrelevant compared to the contributions of u/t-channel for the DM self-interactions for low-velocity systems such as dwarf galaxies. It is evident from the discussion above that the DM-DR interaction cross sections are suppressed by the DM mass. So, for the parameter space of our interest in this work, DM and DR are decoupled very early and cannot significantly change the small-scale structures of the Universe. Before closing the DM section, we would like to mention that we will not study the indirect detection aspects of this model due to the complication with the Sommerfeld enhancement in low-velocity systems.

Electroweak phase transition
The dynamics of the phase transition is determined by the effective potential at the finite temperature (see, e.g., Ref. [111] for a recent review), which can be calculated perturbatively or nonperturbatively on the lattice with dimensional reduction [112][113][114]. While the latter approach provides a gauge independent result and is free of the infrared problem [115], it is computationally expensive and so far has been adopted for only a few models with a simple extended Higgs sector [21,28,44,[116][117][118]. Therefore the perturbative method was predominant in the literature on the analysis of a thermal phase transition. In the standard perturbative approach, the effective potential receives contributions from the tree-level potential, the one-loop Coleman-Weinberg correction and its finite-temperature counterpart, as well as Daisy resummations, which together leads to a gauge dependent result (see, e.g., Refs. [119,120] for a study of the uncertainties with this approach). A gauge independent result nevertheless can still be obtained if only the leading order thermal correction at the high temperature is kept [121]. This also makes an analytical understanding of the otherwise complicated effective potential possible and can better guide the exploration of the phase history. Thus we follow this gauge independent perturbative approach. The finite temperature effective potential can thus be written in the following simplified form where V tree is given in Eq (3.1) and ∆V (1) (T ) is the leading thermal correction given by [122] ∆V (1) where φ i (i = 1, 2, 3) indicates any of the three fields. Here the functions J B and J F have the following high-temperature limit, i.e., for y ≡ m/T 1, Therefore at order y 2 , the thermal corrections reduce to a simpler polynomial form where M S and M V are the field-dependent masses for scalar and dark gauge bosons, which are given in Appendix A. From the finite temperature effective potential, the details of the phase Extrema Type [123] Ref. [123]   Even though those two scenarios have the same zero-temperature potential in Eq. (3.1), the mass parameters evolve differently with temperature as shown in Eqs. (5.5) to (5.10). The parameter space for FOPT in those two scenarios is not the same, though the phase transition pattern should not be qualitatively different. For the rest of this paper, we will focus on the two BMs in Table 2 in the TT scenario as an illustration for the phase transition and GW generation. Given the three possible non-zero VEVs (v h , v 1 , v 2 ), there are eight combinations of possible extrema. Those and their stable conditions are listed in Appendix B and summarized in Table 3. With the desirable features from the extra DR, we require that at T = 0, the stable vacuum be in Type-7: (v h , 0, v 2 ). From the scanning, we found mainly two possible paths of the phase transitions to achieve this pattern two-step: where "⇒" indicates a first-order phase transition and "→" for a continuous transition. 3 The two-step transition as in Eq. (5.11) can yield an electroweak FOPT [20], while the second path in Eq. (5.12) would not lead to an electroweak FOPT and can wash out any previously existing baryon asymmetry. To give a clearer picture of the above transitions, we illustrate the vacuum evolution in detail for the case of BM1 as defined in Table 2. In this case, the phase transition is a two-step process as shown in Eq. (5.11). The evolution of the vacuum and the corresponding potential values in BM1 are shown in the upper panel of Fig. 8. We see that at high temperatures, the stable vacuum is in a symmetric phase of Type-1: (0, 0, 0). At T ≈ 200 GeV, the field Φ 1 develops a VEV and the stable phase becomes of Type-3: (0, v 1 , 0) through a continuous transition, where the order parameter, the VEV v 1 , undergoes a continuous change. As the temperature further decreases, another minimum appears via Φ 2 at (0, 0, v 2 ), which eventually evolves into a minimum of Type-7 (v h , 0, v 2 ) continuously. At T = T c , corresponding to the right of the vertical dashed line in the left panel of Fig. 8, these two types of vacua (0, v 1 , 0) and (0, 0, v 2 ) are degenerate, and are separated by a barrier, characteristic for a FOPT. At T < T c , the initially stable vacuum at (0, v 1 , 0) now becomes metastable while the phase corresponding to (0, 0, v 2 ) becomes energetically preferable, and the Universe becomes supercooled as T decreases. During the coexistence of these two phases, while the probability for the Universe to make a transition from the former to the latter becomes increasingly higher, it remains significantly small during this period. The temperature at which the phase transition happens can be quantified by the temperature when there is about one bubble per Hubble volume, and is called the nucleation temperature T n , corresponding to the left of the vertical dashed line in Fig. 8. As T decreases towards T n , the minimum at (0, 0, v 2 ) evolves into (v h , 0, v 2 ). At T ≈ T n , the transition then proceeds through the formation of bubbles, with the vacuum inside being the more stable one (v h , 0, v 2 ), and that outside the metastable one (0, v 1 , 0). Thus the VEV changes non-continuously. The BM2 has a three-step phase transition shown in Eq. (5.12) and the lower panels of Fig. 8. It is similar to BM1 but different in that it has a prolonged phase at (0, 0, v 2 ) coexisting with the metastable (0, v 1 , 0). The tunneling probability is thus high enough for a FOPT from (0, v 1 , 0) to (0, 0, v 2 ) before the latter evolves into (v h , 0, v 2 ). After this step, the vacuum at (0, 0, v 2 ) makes a further continuous electroweak transition to (v h , 0, v 2 ). Further description of the process is provided in an Appendix C.

Gravitational waves
From studies of the above phase transition and its evolution at different temperatures, one can determine a set of portal parameters that determine the resulting GW signals [124] T n , α e , β/H n , v w , (5.13) where T n , as introduced previously, is the nucleation temperature denoting roughly the time for the onset of phase transition when there is one bubble per Hubble volume; α e is a dimensionless quantity characterizing the energy fraction released from the phase transition in the unit of the total radiation energy density at T n ; β is roughly the inverse time duration of the phase transition determining the peak frequency of the GWs and H n is the Hubble rate H at T n ; v w is the wall velocity. The calculations start with the determination of the tunneling probability per unit time per unit volume given by [125] where S 3 is the three-dimensional Euclidean action corresponding to the critical bubble: with the scalar field minimizing the action and corresponding to the solution of the following equation of motion: which says that there is about one bubble in a Hubble volume. A rough estimation of nucleation temperature T n is usually obtained using the condition S 3 (T n )/T n = 140 [128]. One can further calculate the parameter β where where T * is the GW generation temperature and is approximately equal to the nucleation temperature T n . Similar to the definition of H n , H * is the Hubble rate H at T * . The dimension of β is hertz and it is related to the mean bubble separation at the phase transition(see, e.g., [127,129] for the derivation in Minkowski and FLRW spacetimes), which in turn gives the typical scale for GW production and thus its peak frequency. Moreover, α e is the vacuum energy released from the EWPT normalized by the total radiation energy density where ∆V (T ) = V low (T ) − V high (T ) is the difference between lower and higher phases, and ρ * rad = g * π 2 T 4 /30, g * is the relativistic degrees of freedom at T = T * . For a phase transition in a thermal plasma, as is considered here, the energy released goes in part into the kinetic energy of the plasma, with energy fraction κ v , which sources gravitational waves, and into the heat of the plasma. The flow can also go turbulent, with energy fraction κ turb , which becomes another source for GW production. A fraction of released energy can also go into the gradient of the scalar fields, which however is believed to be of negligible fraction [130] and we will not consider it here.
With these portal parameters, we are ready to calculate the GW energy density spectrum. The GW from a FOPT, as in most cosmic processes, is a stochastic background and can be searched for using the cross correlation method − see recent reviews on theories [124,131,132] and on detection methods [133,134]. It is now generally accepted that there are mainly three sources for GW production during a cosmological FOPT: bubble wall collisions, sound waves, and magnetohydrodynamic (MHD) turbulence. For bubble collisions, the GW is sourced by the stress energy located at the wall and can be understood very well both analytically [135] and numerically [136] under the envelope approximation [137][138][139], where the wall is assumed to be thin and contribution from the overlapped regions is neglected. There has also been recent progress for simulations going beyond the envelope approximation [140][141][142]. However, for a phase transition proceeding in a thermal plasma, it is believed to be of negligible contribution [130]. A significant fraction of the energy released from the phase transition goes to the kinetic energy of the plasma, while the rest heats up the plasma. The kinetic energy of the plasma corresponds to the velocity perturbations of the plasma, which are sound waves in a medium consisting of relativistic particles. This relatively long-living acoustic production of GW is generally accepted to be the dominant one. GW spectrum from this source typically relies on large scale lattice simulations [143][144][145][146]. However, an analytical modeling reproduces the spectra from simulations reasonably well based on the sound shell model [129,147] (see Ref. [127] for the generalization to an expanding Universe), which assumes the plasma velocity field is a linear superposition of the sound shells from all bubbles. The fully ionized fluid can go turbulent for a sufficiently large Reynolds number and corresponds to the third source [143,144]. We will thus include only the contributions from the sound waves and the MHD turbulence, with the present dimensionless GW energy fraction spectrum given by where h ≈ 0.673, the Hubble rate today H 0 in unit of 100 kms −1 Mpc −1 . The sound wave's contribution is [124,148] Ω sw h 2 = 2.65 × 10 −6 H * β κ v α e 1 + α e the spectrum. We note that before the discovery of Υ, a similar suppression factor min(1, τ sw H * ) was adopted [150][151][152] based on a Minkowski derivation of the spectrum [144], which corresponds to the limit of Υ when τ sw H * 1. The lifetime τ sw can be taken as the time scale when the turbulence develops, roughly given by [145,153]: where R * is the mean bubble separation and is related to β through the relation R * = (8π) 1/3 v w /β for an exponential nucleation of the bubbles (see, e.g., Ref. [129] for a derivation in Minkowski spacetime and see Ref. [127] for an analysis in the expanding Universe);Ū f is the root-meansquared fluid velocity and can be determined from the hydrodynamic analysis, with the result U f = 3κ v α/4 [129,148].
The approach above is based on the bag model. As pointed out in a recent work [154], a more robust way to quantify the strenght of phase transition is going beyond the bag model and using the pseudotraceθ to define the strength parameter αθ where the subscript s (b) means evaluating the corresponding variables at the symmetric (broken) phase. The difference inθ in two phases is defined as Dθ ≡θ s (T s ) −θ b (T s ). The speed of sound in the broken phase c s,b is related to the pressure p and energy density e in the broken phase 27) and ω = e + p is the enthalpy density. Beyond the bag model, the GWs power spectrum produced by sound wave can still be obtained from Eq. (5.22) by the replacement α e κ ν α e + 1 → Dθ 4e s κθ. (5.28) In Fig. 9, we show the GWs spectrum within and beyond the bag model and with and without the suppression factor Υ defined in Eq. (5.24). The contributions from MHD turbulence can be modeled as [124] Ω turb h 2 = 3.35 × 10 −4 H * β , (5.29) where κ turb is the energy going to turbulence and f turb is the present day peak frequency: We note that the contribution from MHD is currently the least understood and might witness significant changes in the future. Indeed recent direct numerical simulations show significantly different result [155]. Also the value of κ turb is unknown and we take tentatively κ turb ≈ (5 ∼ 10)%κ v [144] . For both contributions, while in principle the wall velocity v w can be calculated from micro-dynamics of particle interactions with the Higgs condensate, its precise value remains undetermined due to the theoretical uncertainties in the calculations. On the other hand, if baryon asymmetry were to be generated during the phase transition, then a subsonic value is needed. However a supersonic value of v w might still be compatible with EWBG due to the outflowing fluid around the wall [156], as adopted in [45][46][47]157], though a definitive justification of this argument is still missing, which would require a thorough scrutiny of the particle transport near the wall. So we choose tentatively a value v w = 0.9. For the benchmark points in Table 2, the GW spectrum are shown in Fig. 9. To illustrate the suppression effect of Υ, We present the results without considering it by the dashed lines. Some space-based interferometers sensitivities: LISA [158], Taiji [159], TianQin [160], Big Bang Observer (BBO), DECi-hertz Interferometer GW Observatory (DECIGO) and Ultimate-DECIGO [161] are overlaid in Fig. 9. To quantify the detectability of the signals, we define the signal-to-noise ratio (SNR) [124]: where T is the duration of the mission in years. Here we adopt T = 5. The threshold value of SNR for detection is 10 or 50 [124], and thus the BM1 can produce strong GW signal detectable at BBO, though the corresponding values of SNR obtained using the beyond bag model are significantly reduced.

Summary and conclusions
In this work, we extended the Standard Model with a dark SU(2) D gauge sector and a dark scalar sector, beyond the existing literature shown in Table 1. We imposed Z 2 symmetry in the dark scalar sector. We considered two different scalar scenarios under the dark SU(2) D gauge charge, namely, a model with two scalar triplets (TT) and one with a scalar singlet plus a scalar triplet (ST). The dark sector couples to the SM through the Higgs portal − the mixing between the SM Higgs boson and the dark scalars. We worked out the existing constraints on the dark sector model-parameters from the vacuum stability, perturbative unitarity, Higgs physics at the LHC, and the cosmological bounds from CMB measurements and the DM relic abundance and its direct detections. For illustration, we chose two representative benchmark points as shown in Table 2, which satisfy all the constraints, possess the desirable features, and could lead to observable effects. We summarize our novel results as follows: • Higgs boson physics: Via the Higgs portal, the properties of the SM Higgs boson would be modified, including the couplings and an invisible decay. It is particularly interesting to test the potentially large deviation of the Higgs boson triple-self coupling from the SM prediction. Direct searches for the heavy Higgs boson decaying to the SM heavy particles may also be fruitful. We showed those in Figs. 1, 3, and 4.
• DR: Because of the existence of a massless DR associated with the unbroken subgroup U(1) D , it can introduce the velocity-dependent DM self-interaction, which would be desirable to resolve the small-scale structure problems.
• DM: The two stable massive gauge bosons associated with the broken dark gauge group and the pseudo-Goldstone boson can serve as cold DM candidates. The acceptable relic densities were shown in the left panel of Fig. 7. We explored the prospects of their detection in the direct DM searches as shown in the right panel of Fig. 7.
• EWPT: The nontrivial scalar potential has eight types of vacuum pattern for the vacuum structure as shown in Table 3. We have found both the two-step and three-step phase transitions with the cooling of the Universe. Due to the rich vacuum pattern, the scalar sectors can introduce a strong FOPT, as illustrated in Fig. 8 for the benchmark points BM1 with a successful EW FOPT, and BM2 with a FOPT in the dark sector.
• GW: Our benchmark GW spectra are shown in Fig. 9. We found that the two-step EWPT in our BM1 can produce strong GW signals and can be detectable using the future spacebased interferometers BBO, while the GW signal for BM2 may be difficult to observe at BBO due to the rather low signal-to-noise ratio. We summarize our results on the m h 2 -sin θ plane in Fig. 10, fixing the other parameters according to our BM1 (left panel) and BM2 (right panel). The orange shaded regions are allowed by the DM direct detections. Outside the cyan shaded regions, DM would over-close our Universe. The black points are the viable FOPT points which can enable GW production. The gray solid lines show the predicted deviation of the SM triple Higgs coupling. Our BM1 and BM2 points sit in the red-cross and blue-star symbols, respectively.
In conclusion, given the outstanding puzzles we are facing now such as the identity of the DM and the nature of the EWPT, it is prudent to consider the possibility of a dark sector uncharged under the SM interactions. We demonstrated with a well-motivated example of a dark SU(2) D sector, that rich physics may exist that is potentially observable with the current and future measurements at colliders, DM experiments, and GW interferometers.

Appendices A Field-dependent mass
The field-dependent masses for the scalar degrees of freedom with n S = 1 are Similarly the field-dependent masses for the vector degrees of freedom with nW = n Z = 3 and n W = 6 are m 2W 1 =g 2 (ϕ 2 1 + ϕ 2 2 ), m 2W

B Stable conditions for all the minima
The stable conditions for the extrema in Table 3 are Those cases are summarized in Table 3. For Type-8 scenario, we refer to Ref. [123] due to the complicity and irrelevance.

C Further description for the phase transition process
The effective potential is a polynomial of the fields (h 0 , ω 3 , ϕ 2 ) up to quartic terms by renormalizibility. The coefficients of the quartic terms are required to be positive as the potential is bounded from below. For the quadratic terms, they can be generically put into the form with φ i denoting one of the three fields. For D i > 0 and T > T i , this term remains positive and enforces a minimum at φ i = 0, which is in a symmetric phase. As T decreases below T i , the minimum at φ i = 0 will roll away from the origin and takes a non-zero value, corresponding to a continuous phase transition. This is indeed what happens for the continuous transitions in Eq. (5.11) and Eq. (5.12), where the potential minimum corresponds to a non-zero field value at some temperature. The same story can happen to any of the three fields. If the parameters are such that two minima coexist across a time duration, then a first order phase transition can happen when the universe tunnels from one minimum to another, characteristic for a firstorder phase transition (FOPT), as shown in these two benchmarks in the text. This analytical understanding can provide a way of identifying the parameter space giving a first order phase transition, as demonstrated in [50,76]. In practice, however, there is a challenge in this procedure. Whether or not a transition takes place between two coexisting minima depends on the tunneling probability and it is sensitive to the potential shape such as the height of the barrier separating them and the potential difference at the two minima, which however is difficult to understand analytically (see [46,162] for relevant analyses and discussions). This presents an uncertainty for the presence of a FOPT even if we perceive the coexistence of two minima at the same time.
As such, some numerical techniques, such as scanning over a large parameter space, may be unavoidable, as we did in our analyses. There are also subtleties in classifying second-order/higher-order phase transitions and a smooth cross-over. A proper classification could be specified by a dimensionless susceptibility, see, e.g. Ref. [118].